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ABSTRACT 

The  nutational  stability  of  a  dual-spin,  quasi-rigid,  axisymmctric  spacecraft  containing 
a  driven  rotor  is  analyzed.  The  purpose  is  to  examine  a  revised  energy-sink  stability  theory 
that  properly  accounts  for  the  energy  contribution  of  the  motor.  An  inconsistency  in  the 
development  disproves  the  existing  energy-sink  theory's  assumption  that  the  motor  of  the 
system  contributes  exactiy  enough  energy  to  offset  the  frictional  losses  between  the  rotor 
and  the  platform.  Using  the  concept  of  core  energy,  the  revised  stability  criteria  for  a  dual- 
spin,  quasi-rigid,  axisymmetric  spacecraft  containing  a  driven  rotor  is  derived.  An 
expression  for  nutation  angle  as  a  function  of  core  energy  over  time  is  then  determined. 
Numerical  simulations  are  used  to  verify  the  revised  energy-sink  stability  theory.  The  dual- 
spin,  quasi-rigid,  axisymmetric  system  presented  by  D.  L.  Mingori  was  chosen  for  the 
simulation.  Equations  for  angular  momentum  and  total  energy  were  necessary  to  validate 
the  numerical  simulation  and  confirm  aspects  of  the  revised  energy-sink  stability  theory. 
These  equations  are  derived  from  the  first  principles  of  dynamics  and  are  included  in  the 
analysis.  An  explicit  relationship  for  core  energy  as  a  function  of  time  does  not  exist. 
Various  models  postulating  core  energy  are  presented  and  analyzed.  The  numerical 
simulations  of  the  computed  nutation  angles  as  a  function  of  the  postulated  core  energy 
compare  well  with  the  actual  nutation  angles  of  the  system  to  confirm  the  revised  energy- 
sink  stability  criteria. 
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I.  INTRODUCTION 

A  chronological  review  of  early  spacecraft  types,  and  the  stability  criteria  developed 
for  them,  provides  a  sufficient  background  for  the  fundamental  concepts  of  this  thesis.  The 
revised  energy- sink  stability  criterion  is  then  presented,  and  an  equation  for  nutation  angle 
as  a  function  of  energy  dissipation  over  time  is  developed. 

A.    SINGLE  SPIN  SATELLITES 
1.  Equations  of  Motion 

The  earliest  satellites  took  advantage  of  the  faa  that  stability  could  be  achieved  via 
the  'gyroscopic  stiffness'  of  a  spinning  body.    A  preliminary  dynamic  model  for  the 
satellite  could  be  achieved  with  Euler's  equations  of  motion  for  a  rigid  body.   Eu' 
moment  equation  can  be  written  as 

M  =^ih  =^^  +  ^a)^xh  (1, 

dt         dt 

In  component  form  this  becomes 

Af  1  =  hi  +  (02  hs  -  ci>2  /z2 

Af2  =  ^2  +  0)3  h\  -  0)1  /l3  (2) 

M3   =   h3  +  CDi  /l2  -  0)2  hi 

Simplification  of  the  model  is  accomplished  by  assuming  that  it  is  axisymmetric,  that  the 
body-fixed  axes  coincide  with  the  principal  axes  (correct  for  simple  spacecraft  with 
/i  =  I2),  and  that  the  body  is  in  a  torque-free  environment  (a  valid  approximation  used 
throughout  this  thesis).  The  spacecraft  is  then  represented  by  Figure  (1)  and  the  equations 


of  motion  are 


0  =  /i  ODi  +  {I3  -  Ii)  0)2  Q>3 

0  =  /2G)2  +  (/l-/3)fl>3^1 
0  =  /3  G)3 


(3) 


l>3 


Axisymmetric  Rigid  Body 
Figure  (1) 


The  angular  velocity  vector,  angular  momentum  vector,  and  the  kinetic  energy  may  be 
expressed  respectively  as 


CD     =   G>i  bi  +  0)2  b2  +  0)3  bs 

h  =  /i  Q>i  bi  + 12  fl>2  b2  + 13  C03  b3 
E  =  ^(licof  +  hcoi  +  hco^) 


(4) 
(5) 

(6) 


A  simplification  of  the  notation  can  be  made.  Let  /i  =  /2  =  A  and  h^Is  •  From  the  third 
line  of  Equation  (3)  it  can  be  seen  that  the  angular  velocity  about  the  spin  axis  b^  is 
constant,  therefore  03  =  0)5  .  The  transverse  angular  velocity  components  interchange 
between  the  bi  and  the  b2  axes  but  the  magnitude  is  constant,  so  one  can  let 


(Of  =  fi)i  bi  +  G)2  b2.  Therefore, 

0)     =  COr  +  CDj  b3  (7) 

h  =  It(ih  +  Is(0s^3  (8) 

/|2  =  //  0,^2  ^  ;^2  ^^2  (9) 

Because  the  motion  is  torque-free,  |h|  is  constant  and  h  is  fixed  in  space.  Because  there  are 
no  energy  sources  or  sinks,  E  is  also  constant.  Finally,  the  above  conditions  result  in 

I  CO  I  being  a  constant. 

An  expression  for  the  nutation  angle  may  now  be  developed.  The  orientation  of 

the  body  axes  with  re;       t  to  an  inertial  reference  frame  is  desired  to  provide  a  measure  of 

the  body's  dynamic  behavior.  Because  the  angular  momentum  vector  is  fixed  in  space,  the 

nutation  angle  is  defined  as  the  angle  between  the  body-fixed  axis  about  which  spin  is 

desired  and  the  angular  momentum  vector,  and  can  be  expressed  in  one  of  the  following 

forms 

e  =  tan-^  f^^^Vr/r^l  =  tan-^  f^i  ^i  »>i  ^2  a>2  bjj  ^  ^^_,  j^j 

\         |/I3b3|         /  I  h0)3b3\  j  UsCOsI 

In  these  equations  the  first  and  second  terms  correspond  to  the  general  case  with  spin  about 
the  bs  axis,  and  the  third  term  uses  the  simplified  notation  to  describe  an  axisymmetric 
body.  In  the  special  case  of  spin  about  only  one  principal  axis  of  an  axisymmetric  body, 
the  angular  momentum  vector  and  the  angular  velocity  vector  will  lie  on  the  spin  axis.  With 
spin  components  on  two  or  more  principal  axes,  the  three  vectors  will  not  be  coincident, 
although  they  still  will  lie  on  the  same  plane. 


2.  Single  Spin  Satellite  Stability 

A  torque-free,  axisymmetric  rigid  body  with  the  body  axes  coinciding  with  the 
principle  moments  of  inertia  will  be  stable  about  the  axis  of  either  the  maximum  moment  of 
inertia  or  the  minimum  moment  of  inertia.  To  prove  this,  one  begins  with  an  arbitrary  rigid 
body.  The  body  is  given  the  initial  condition  of  steady  angular  velocity,  coq,  about  a 
principal  axis,  and  is  then  perturbed  slightiy.  It  is  assumed  that  the  angular  velocities  about 
the  other  axes  are  small,  and  are  approximately  the  same  order  of  magnitude  as  the 
perturbation  (coi  =  6)2  =  £}•  The  system  will  be  considered  stable  if  the  perturbation  does 
not  increase  over  time.     Given  an  initial  angular  velocity  with  a  perturbation, 

CO   =  CO  I  hi  +  C02b2+{o)o  +  e)b2,  and  given  arbitrary  inertias  /i,   I2,    h^   Euler's 
equations  of  motion  can  be  written  as 

0  =  /i  6)1  +  [h  - 12)  0)2  {coq  +  £) 

0  =  /2  a>2  +  (/i  -  h)  (Oi  (0)0  +  e)  (14) 

0  =  /3  0)3  +  (/2  -  /i)  0)i  0)2 

TTie  equations  are  linearized  by  neglecting  terms  of  magnitude  £^.  Rewriting  the  equations 
by  eliminating  the  terms  ecoi,  £  0)2,  and  0)i  0)2  results  in 

6)1  =  ^^^^  0)2  0)0  (15) 

(/3-/1) 

0)2  =      \         0)1  0)0 

h  (16) 

0)3    =    6)Q  +  £   =    £    =    0  /jyx 

From  Equation  (17),  one  can  conclude  that  £  =  constant.  By  differentiating  Equation  (15), 
and  using  Equation  (16)  to  eliminate  0)2,  one  gets 

Similarly,  from  differentiating  Equation  (16),  and  using  Equation  (15)  to  eliminate  0)i,  one 


amvesat 


^^jiJlSzlMlldAtoi^c^^O  (19) 

These  equations  are  identified  as  second  order,  linear,  ordinary  differential  equations  with 
constant  coefficients.  The  general  solution  for  these  differential  equations  are 

0)2  =  K3tirf  +  K4t-irt  (21) 

where 


■V 


/l  /2 


If  >  is  imaginary,  coi  and  CO2  will  increase  without  bound  over  time,  and  the  motion  will  be 
unstable.  Stability  is  achieved  if  y  is  real.  The  first  case  occurs  when  the  maximum 
moment  of  inertia  is  about  the  spin  (bs)  axis;  then  /s  >  /i  ,h>  h*  and  {Iz  -  h)  {h  -  ^1 )  >  0- 
In  this  case  the  inertia  ratio  Is  I  It,  defined  as  the  inertia  about  the  spin  axis  over  the  inertia 
about  a  transverse  axis,  is  greater  than  one.  The  second  case  occurs  when  the  minimum 
moment  of  inertia  is  about  the  spin  (bs)  axis;  then  h<I\  ,  /a  <  h,  and  (/s  - 12)  (/s  -  /i)  >  0 
as  well.  Here  the  inertia  ratio  is  less  than  one.  In  both  of  the  above  cases  y  is  real  and  the 
motion  is  stable.  However,  if  /s  is  the  intermediate  moment  of  inertia  about  the  spin  (ba) 
axis,  then  [I^  - 12)  {I3  -  /i)  <  0,  y  is  imaginary,  and  the  motion  is  unstable. 

The  previous  model  cannot  be  applied  to  a  satellite  since  the  assumption  of  a  rigid 
body  cannot  be  extended  to  the  spacecraft.  Structural  elasticity,  liquid  propellant  slosh, 
etc.,  cause  energy  dissipation  in  an  actual  spacecraft.  This  spacecraft  can  be  generalized  by 
a  quasi-rigid  body  with  an  unspecified  energy  damper  mechanism.  A  priori,  one  can 
conclude  that  energy  in  the  above  system  will  dissipate  until  the  minimum  energy  state  is 
reached.  The  kinetic  energy  of  the  system  with  spin  about  the  principal  axis  with  maximum 
moment  of  inertia  and  spin  about  the  principal  axis  with  minimum  moment  of  inertia  can  be 


written  respectively  as 


E^\-b—  /;^aboutb3 

^  'max 

(23) 
£  =  i-^  /min  about  b3 


■min 


The  angular  momentum  is  constant  in  the  torque-free  case.  Thus  the  minimum  kinetic 
energy  state  occurs  when  the  rotation  is  about  the  axis  of  the  maximum  moment  of  inertia. 
Thereforc,  a  quasi-rigid  body  is  stable  only  when  it  is  spinning  about  its  major  axis,  with  a 
corresponding  inertia  ratio  that  is  greater  than  one. 

A  relationship  can  be  established  between  the  time  rate  of  change  of  the  nutation 
angle  and  the  energy  dissipation  of  a  quasi-rigid,  axisymmetric  body.  One  must  assume 
that  the  angular  momentum  and  moments  of  inertia  of  the  quasi-rigid  body  do  not  change 
appreciably  from  a  comparable  rigid  body.  For  the  generalized  model,  with  arbitrary 
inertias  I\  =  I2,  and  /s.  Equations  (5)  and  (6)  are  substituted  into  Equation  (12)  to  obtain 

sin2  e   =  (_A_){2^^Z^  (24) 

The  time  rate  of  change  of  the  nutation  angle  is  determined  by  taking  the  derivative  of  the 
above  equation 

e  =  —j—r      ^'^^\  Ed  total  (25) 

where  the  only  rate  of  change  of  energy  E  is  attributed  to  the  damping  mechanism  and  is 
written  as  Ep  total  to  emphasize  this  point.  Because  Eq  total  is  negative,  the  nutation  angle 
will  decrease  only  if /a  is  greater  than  /i.  This  reaffirms  the  previous  conclusion  that  a 
quasi-rigid  body  is  spin  stabilized  only  about  the  axis  of  maximum  moment  of  inertia.  The 
foregoing  development  is  referred  to  as  the  energy-sink  method- 


B.    DUAL  SPIN  SATELLITES 

The  logical  progression  from  the  single  spin  satellite  was  to  incorporate  a  de-spun 
platform.  This  permitted  the  replacement  of  the  low  gain  omnidirectional  antenna  with 
directional  antennas  for  communication  satellites,  and  a  more  capable  spacecraft  for 
scientific  observation.  A  simple  control  system  about  the  bs  axis  would  maintain  the 
platform  rotating  at  a  constant  relative  rate  with  respect  to  the  rotor  (  and  would  usually 
have  the  platform  rotate  at  the  earth's  rotational  rate).  Initially,  the  platforms  were 
sufficientiy  small,  and  the  overall  dimensions  of  the  satellite  were  such  that  the  inertia  ratio 
would  be  greater  than  one.  For  this  type  of  satellite,  the  previously  developed  theory 
proved  adequate.  However,  as  satellites  continued  to  grow  in  size,  the  launch  vehicle 
shroud  diameter  became  a  constraint.  In  order  to  provide  the  size  spacecraft  needed  to 
satisfy  mission  requirements  and  still  fit  within  the  shroud,  a  spacecraft  with  an  inertia  ratio 
of  less  than  one  (I $  total  I  h  total  <  1)  would  need  to  be  built.  From  the  previously 
developed  theory,  a  spacecraft  with  an  inertia  ratio  of  less  than  one  was  believed  to  be 
inherently  unstable.  It  was  not  until  the  development  of  the  energy-sink  theory  for  dual- 
spin,  quasi-rigid,  axisymmetric  spacecraft  containing  a  driven  rotor,  developed 
simultaneously  by  V.  D.  Landon  (unpublished  work)  and  A.  J.  lorillo  [Ref.  1],  that  a 
spacecraft  with  an  inertia  ratio  less  than  one  was  considered  feasible.  Several  rigorous 
stability  analyses  using  the  equations  of  motion  for  specific  dual  spinners  have  been 
performed  by  P.  W.  Likins  [Ref.  2],  D.  L.  Mingori  [Ref.  3],  and  others,  to  validate  the 
energy-sink  theory.  The  difficulty  of  a  rigorous  analysis  is  in  accurately  modeling  all  the 
forms  of  energy  dissipation.  A  more  general  and  practical  approach  was  required  to 
determine  stability,  and  the  energy-sink  theory  proved  suitable.  The  development  of  this 
theory  is  as  follows. 


The  spacecraft,  shown  in  Figure  (2),  is  assumed  to  fulfill  the  following  conditions: 

•  Both  the  rotor  and  the  platform  arc  axisymmetric 

•  Both  the  rotor  and  the  platforai  arc  quasi-rigid 

•  The  damping  mechanisms  do  not  significantly  alter  the  energy  value,  although 

the  mechanisms  will  contribute  to  the  energy  rate 

•  No  external  torques  arc  applied 

•  The  only  rclative  motion  is  spin  about  the  bs  axis 

•  The  motor,  driven  by  the  control  system,  inputs  just  enough  energy  to  exactly 

offset  the  shaft  frictional  losses,  maintaining  a  constant  relative  angular 
velocity  between  the  rotor  and  the  platform 


»>3.b'3.b"3  i 


Dual-Spin  Quasi-Rigid  Axisymmetric  Spacecraft 
Figure  (2) 
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The  magnitude  of  the  angular  momentum  and  the  kinetic  energy  of  the  dual-spin  system  can 
be  expressed  respectively  by  the  following  equations 

h'^   =  J^total  Oi?  +  [h  COs  +  1/  COsf   =   It^total  ^/^  +  (^5  CT  +  Is  total  CO/f  (26) 

£^  =  |//  total(0i^+^IsC0s^+^I/C0/^=-^It  totalOt^+^Is  totalCO/^ +jIsO'^+IsC0/o       (27) 

where  G  is  the  relative  rotation  rate  of  the  rotor  with  respect  to  the  platform.  Although  the 
equations  actually  represent  a  rigid  body  system,  they  are  also  applicable  to  the  quasi-rigid 
system  because  of  the  above  assumptions.  If  the  damping  mechanisms  do  make  significant 
contributions  to  the  energy  of  the  system,  then  Equations  (26)  and  (27)  do  not  hold,  and 
the  energy-sink  criterion  will  not  apply.  Additionally,  the  potential  energy  of  the  system 
(for  example  the  energy  stored  in  the  spring  of  a  mass-spring-dashpot  damper)  will  not 
make  a  significant  contribution  to  the  total  energy  of  the  system.  Therefore,  for  the 
remainder  of  the  thesis,  it  is  assumed  that  the  kinetic  energy  of  the  system  is  effectively  the 
total  energy  of  the  system  (E  =  Etotaih  Because  the  system  experiences  no  exiemal  torques, 
angular  momentum  is  conserved.  Because  the  motor  contributes  no  energy  to  the  system, 
the  time  rate  of  change  of  the  kinetic  energy  E  is  represented  by  only  the  quasi-rigid  energy 
dissipation  Ep  totaU  and  is  negative.  Differentiating  the  above  two  equations  with  respect  to 
time,  one  obtains 

0   =   tuotal  ^r  ^/  +  [h  COs  +  1/  CO/)  (h  0)s  +  V  d)/)  (28) 

E  =  Ed  total  =  A  total  COt  fi)/  +  Is  0)s  4  +  1/  0)/  6)/  (29) 

Eliminating  the  common  term  cOt  d>t  by  combining  the  above  equations 


^  .  -\Is  COs  +  //  (Oslils  (Os  +  Is'  CDs')       ,  ,,,'   ,  ,^«v 

E  =  Ed  total  =  -^^^—^ —     /  — '—^  + Is  COs  COs +  1/0)/ CO/        (30) 

't  total 

One  may  now  define  the  inertial  nutation  frequency,  Xo,  the  rotor  nutation  frequency. 


A,  and  the  platform  nutation  frequency.  A',  as  follows 


'  ,^ ' 


.  IgCOs  +  Is   CDs  ,^,, 

Ao  =  ^-^ — - — -  (31) 

't  total 

^         o         ,,  (Is-Ittotal)(Os-^Is   (Os  .-^. 

A  =  A,Q-(Os  = (32) 

h  total 

A    =  Ao  -  ^5    =  7 (33) 


The  nutation  angle  for  a  dual- spin  system  is  defined  as 

^  Jh>b3^  i(^3+^\  ,f/5  6?5-H//Q)/ 


(34) 


By  imp>osing  the  condition 

A^  >  0  (35) 

the  analysis  of  nutational  motion  is  restricted  to  the  following  region  without  any  loss  of 
generality 

0  <  0  <  I  (36) 

Incorporating  the  nutation  frequency  terms,  the  equation  for  energy  dissipation  is  written  as 

E  =  Ed  total  =  Ed  +  Ed'  =  -  /,  A  G)^  -  Is'  A'  6)s'  (37) 

Recalling  the  assumption  that  the  rotor  and  platform  arc  uncoupled  about  the  bs  axis,  we 
may  incorporate  Equation  (37)  into  the  reaction  torques  which  tend  to  change  the  angular 
rates 

^  (38) 

X' 
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Combining  Equation  (30),  (37),  and  (38),  one  arrives  at  the  transverse  rate  equation 

It  total  COt  6)t  =  -^  — +  — I 


(39) 


Differentiating  Equation  (12)  and  substituting  Equation  (39)  into  it,  the  time  rate  of  change 
of  nutation  angle  as  a  function  of  energy  dissipation  rates  is 


Q  =  2  It  total  Ac  JEp  ^  Ep'^ 
sm[2e)h?-\X       V 


(40) 


The  energy-sink  equation  for  a  single  spin-stabilized  body  is  obtained  by  letting  A  and  A' 
become  Ao  and  Ep  +  Ep  become  Ep  total  for  a  single  body.   The  definition  of  stability 

requires  the  nutation  angle  Q  to  remain  constant  or  decrease  as  a  function  of  time,  so  that  Q 
is  zero  or  negative.  Because  the  factors  outside  the  parenthesis  on  the  right  hand  side  of 
Equation  (40)  are  positive,  the  stability  criterion  for  a  dual-spin,  quasi-rigid,  axisymmetric 
system  becomes 

(41) 


One  of  the  following  cases  will  guarantee  stability 

1)     A>0   and   A'>0 


2)     A>0,  A'<0,  and 


3)     A<0,  A'>0,  and 


A  specific  example  would  be  the  model  of  a  typical  communications  satellite,  a  prolate 
dual-spinner  possessing  an  inertia  ratio  of  less  than  one.  In  general,  the  rotor  nutation 


Ep 

< 

Ep 

X' 

A 

Ep' 

> 

Ep 

r 

A 
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frequency,  expressed  as 

_  Is  COs  +  1/  CO/      ^^     _  {Is  -  It  total)  COs  +  //  CO/ 

A   —   J COs   —  7  V*^) 

h  total  *t  total 

would  be  negative  because  {Is  -  It  total)  is  negative  and  cOs  »  cOs'  if  cOs'  is  rotating  at  the 
earth's  rotation  rate.  Thus,  energy  damping  in  the  rotor,  Ep,  would  be  destabilizing  while 
energy  damping  in  the  platform,  £"£>',  would  be  stabilizing.  It  is  from  this  result  that 
satellites  will  have  a  damping  mechanism  placed  on  the  platform  to  improve  nutational 
stability.  Such  a  damper  is  called  a  nutation  damper. 
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n.  PROBLEM  DEFINITION 

A.    OVERVIEW 

The  existing  energy-sink  theory  relies  on  several  assumptions,  perhaps  the  most 
important  relating  to  the  driven  rotor.  As  previously  stated,  it  has  been  assumed  that  the 
motor,  driven  by  the  control  system,  inputs  just  enough  energy  to  exactiy  offset  the  shaft 
frictional  losses,  maintaining  a  constant  relative  angular  velocity  between  the  rotor  and  the 
platform.  In  actual  systems,  contrary  to  this  assumption,  the  motor  may  add  or  remove 
energy  from  the  system,  depending  on  the  dynamics  of  the  spacecraft.  Consequently,  a 
revised  energy  sink  stability  theory,  properly  accounting  for  the  energy  contribution  of  the 
motor,  is  derived.  The  revised  theory,  based  on  the  concept  of  core  energy,  will  remain 
consistent  with  the  existing  energy-sink  stability  criterion.  Continuing,  an  equation  for 
nutation  angle  over  time,  as  a  function  of  core  energy,  is  developed.  Given  a  postulated 
energy  dissipation  function  modeling  the  nutation  dampers,  structural  elasticity,  fuel  slosh, 
etc.,  one  can  accurately  predict  the  nutation  angle  behavior.  Numerical  simulation  of  D.  L. 
Mingori's  dual-spin,  quasi-rigid,  axisymmetric  system  containing  a  driven  rotor  [Ref.  3]  is 
used  to  validate  the  revised  energy-sink  stability  theory.  The  predicted  nutation  angle, 
based  on  this  revised  energy-sink  theory,  and  the  postulated  energy  dissipation  function,  is 
compared  to  the  exact  nutation  angle  of  the  Mingori  system.  By  using  a  suitable  postulated 
energy  dissipation  function,  one  can  achieve  exceUent  agreement  between  the  predicted  and 
the  exact  nutation  angle.  I.  Michael  Ross  [Ref.  4]  performed  this  analysis  on  a  dual  spin 
system  with  a  damper  on  the  platform  only.  The  remainder  of  this  thesis  will  use  the  same 
analysis,  but  on  the  Mingori  system  with  dampers  on  both  the  rotor  and  the  platform. 
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B.  INCONSISTENCY  OF  THE  ENERGY  SINK  THEORY 

There  exists  a  contradiction  between  the  existing  energy-sink  theory  and  the  nutation 
angle  derived  from  it.  This  will  provide  the  motivation  for  developing  a  revised  energy- 
sink  theory  and  an  alternative  equation  for  nutation  angle  over  time.  From  before,  the 
existing  energy-sink  stability  criterion  can  be  expressed  as 

(41) 


and  the  nutation  angle  for  the  dual-spin  system  was  defined  as 


(34) 


As  previously  stated,  for  a  prolate  dual-spinner  {Is  total  I  h  total  <^\  energy  dissipation  in 
the  platform  is  stabilizing  and  energy  dissipation  in  the  rotor  is  destabilizing.  The  angular 
velocity  of  the  platform  about  the  spin  axis,  (Os\  can  be  expressed  in  terms  of  the  angular 
momentum  and  the  kinetic  energy  of  the  system.  Combining  Equation  (26)  and  Equation 
(27),  one  arrives  at 


(O;   =    ^5  g  ^   //  /^  g  j2  h^-2I,totalE  +  Is  C^  {It  total  -  h)  ^^2) 

h  total        V  Ms  total'  I s  total  Vt  total  ~  Is  total) 

Substituting  this  expression  into  Equation  (34)  results  in 

d   =  cos-i(±Vg)  (43) 

where  Q  is  represented  as 

Q  =  (2  £  -  /^  g ')  Is  total  +  I]  g^  h  total  (44) 

l2  / 1  _  IsjotgA  h  total  -  h  total 

^        It  total' 

Initial  conditions  at  r  =  0  will  determine  the  correct  sign,  with  continuity  considerations 
maintaining  the  sign  for  all  of  r  >  0.  Additionally,  no  external  torques  are  applied  to  this 
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system,  resulting  in  |h|  being  constant  Differentiation  of  Equation  (43)  results  in 

e  =  -    .       /     _,  f  (±Ve)  =  ±-S^-L^{ziojodluo!^  (45) 

V  'r  total' 
From  the  definition  of  nutation  angle.  Equation  (34),  and  the  condition  imposed  on  it. 
Equation  (36),  the  positive  sign  must  be  chosen  in  Equations  (43)  and  (45).  An  important 
observation  is  made  at  this  time.  Choosing  the  positive  sign  will  result  in  a  positive  rate  of 
change  in  the  nutation  angle,  indicating  an  unstable  condition.  The  relative  rotor  spin  rate  is 
an  independent  variable,  and  is  arbitrarily  selected  here  as  a  constant  value  over  time. 
Therefore,  energy  dissipation  in  a  prolate  dual-spin  spacecraft  will  make  the  nutation  angle 
increase,  regardless  of  whether  the  dissipation  is  in  the  platform  or  in  the  rotor.  This  is  not 
consistent  with  the  stability  criterion  of  Equation  (41).  Thus,  the  existing  energy-sink 
stability  criterion  contradicts  itself. 

C.    CORE  ENERGY  AND  ENERGY  DISSIPATION 

The  existing  energy-sink  stability  criterion  does  not  properly  account  for  any  energy 
that  may  be  provided  by,  or  absorbed  by,  the  motor.  To  accurately  represent  the  system, 
the  total  rate  of  change  of  energy  must  be  written  as 

E  =  Ed  total +W  (46) 

where  W  is  the  rate  of  work  due  to  the  motor,  and  may  be  either  positive  or  negative  and 
the  rate  of  change  of  energy  due  to  dissipative  elements  can  occur  in  either  the  platform  or 
the  rotor.  Recall  that  the  kinetic  energy  of  the  dual-spin  system  was  expressed  in  Equation 
(27).  If  the  work  due  to  the  motor  torque  as  a  function  of  time  is  written  in  analytical  form, 
then  the  time  rate  of  change  of  the  energy  of  the  system  due  to  all  dissipative  elements, 

Ed  totaiy  can  then  be  expressed  solely  in  terms  of  the  quasi-rigid  parameters  of  the  dual-spin 
system.   With  this  expression,  the  condition  that  Ed  total  ^  0  will  result  in  the  required 
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stability  criterion.  The  difficulty  arises  in  that  to  get  the  work  due  to  the  motor  torque  W 
(or  W),  one  needs  to  know  the  exact  dynamics  of  the  dissipative  mechanism.  In  the 
development  of  the  modified  energy-sink  stability.  Equation  (46)  will  be  used  to  determine 
the  expression  for  Ep  totah  and  ultimately  derive  the  revised  stability  criterion  and  nutation 
angle  equation. 

Additionally,  the  existing  energy-sink  theory  can  be  shown  to  be  incorrect  for  both 
the  case  of  total  energy  decreasing  and  for  the  case  of  total  energy  increasing.  For 
example,  if  the  rotor  is  rigid  and  total  energy  decreases.  Equation  (41)  predicts  that  the 
system  will  be  stable,  but  Equation  (45)  predicts  that  it  will  be  unstable.  Allowing  the  total 
energy  to  increase  would  reverse  the  conditions,  but  still  show  a  contradiction  between 
Equation  (41)  and  Equation  (45). 

A  modification  of  the  energy-sink  stability  theory  and  the  associated  expression  for 
nutation  angle  is  now  presented.  The  development  of  the  theory  is  from  I.  Michael  Ross' 
unpublished  notes.  The  basis  of  the  new  theory  is  centered  on  the  core  energy  of  the 
system.  As  defined  by  Hubert  [Ref.  5] 

Core  energy  is  the  total  energy  of  the  spacecraft  minus  the  portion  of  the  rotor 
energy  that  is  due  to  the  relative  rotation  between  the  rotor  and  the  platform.  It 
is  assumed  that  the  mass,  inertia,  and  motion  of  the  damping  device  are 
sufficiently  small  that  its  energy  is  negligible  relative  to  the  spacecraft  core 
energy.  The  damper  will  be  treated  as  an  undefined  'energy  sink'  for  the 
purposes  of  the  energy  sink  analysis. 

From  the  above  statement  one  can  define  the  core  body  as  that  body  whose  inertial 
dynamics  are  selected  for  analysis.  Hubert  defines  the  platform  as  the  core  body.  The  core 
energy  is  simply  the  rotational  kinetic  energy  of  a  fictitious  rigid  body  that  possesses  the 
inertia  properties  of  the  entire  dual-spinner  but  moves  in  inertial  space  exactly  like  the  core 
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body.  For  a  dual-spin  spacecraft  the  platform  core  energy  is  defined  as 

Ec    =  ^  /l  total  0>l+^h  total  0)i  +  ^ h  total  0>3^   =  ^  h  total  (^}  +  ^  ^5  total  Ois         (47) 

The  center  expression  is  for  an  arbitrary  dual-spin  spacecraft  with  the  spin  axis  about  the  bs 
axis,  and  the  right  expression  is  for  a  dual-spin,  axisymmetric  spacecraft  with  the 
simplified  notation.  Extending  this  concept  to  the  rotor,  the  rotor  core  energy  is  defined  as 

Ec  =  ^htotal(Oi   +^l2  total  Coi  +  ^h  total  0)3^   =  ^  It  total  COf^  +  ^ Is  total  COg'^  (48) 

Necessary  to  the  development  of  the  modified  theory  is  what  will  be  termed  the 
Separation  Axiom.  This  is  when  analysis  is  first  performed  with  the  rotor  considered  rigid 
and  the  platform  quasi-rigid.  Euler's  equations  are  written  for  the  rotor,  and  through 
manipulation,  an  equation  is  derived  relating  the  torque  on  the  quasi-rigid  platform  solely  in 
terms  of  platform  variables.  Then  the  platform  is  considered  rigid  and  the  rotor  quasi-rigid. 
An  equation  is  derived  relating  the  torque  on  the  quasi-rigid  rotor  solely  in  terms  of  rotor 
variables.  These  two  separate  equations  are  then  combined  and  applied  to  a  system  in 
which  the  rotor  and  the  platform  may  both  be  quasi-rigid. 

The  case  of  a  rigid  rotor  with  a  quasi-rigid  platform  is  first  analyzed.  Because  the 
rotor  is  rigid,  the  torque  applied  by  the  motor  to  the  rotor  is  the  net  torque  on  the  rotor  and 
is  determined  from  Euler's  moment  equations.  For  the  case  of  the  axisymmetric  rotor, 

Tr  =  Tr/m  =  Is  0)s  (49) 

One  can  observe  that  Tp/M  =  -  Tr/m  ,  but  Tp  ^  Tr  since  the  damping  mechanism 
contributes  additional  torques  to  the  quasi-rigid  platform.  The  rate  of  work  needed  by  the 
motor  torque  to  maintain  constant  relative  motion  between  the  rotor  and  the  platform  is 

W  =  TrC  =^  IsGQ)s  =  Is  c(o)/  +  or)  (50) 

By  substituting  the  platform  core  energy.  Equation  (47),  into  the  kinetic  energy  expression. 
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Equation  (27),  kinetic  energy  may  be  expressed  as 


(51) 


The  above  equation  is  differentiated  to  arrive  at  the  time  rate  of  change  of  the  kinetic  energy 
of  the  dual-spin  system 

E  =  Ec  +Is(yc^  +  ^s (era)/  +  cOs'  cj  =  Ec  +  Is  ^(cOj'  +  o)  +  /,  cOs'  6  (52) 
Substituting  into  this  equation  the  rate  of  work  done  by  the  motor  torque,  Equation  (50), 
one  gets 


E  =  Ec+W  +  Is  CO/  6 
Comparing  this  to  Equation  (46),  it  can  be  seen  that 


(53) 


Ed  total  =  Ec  +  h  0)/  6  (54) 

Taking  the  derivative  of  Equation  (47)  to  get  the  time  rate  of  change  of  the  platform  core 
energy 

Ec   =  ft  total  COtQ)t  +  Is  total  COs'  O^s'  (55) 

and  then  substituting  this  into  Equation  (54),  one  arrives  at  the  total  energy  dissipation  of  a 
rigid  rotor,  quasi-rigid  platform  system 


Ed  total  =  It  total  (Ot6)t  +  Is  total  0>s'  C^s'  +  h  G>/  O" 

=  It  total  COtQ)t  +  Is  CO/  6)/  +  //  CO/  6)/  +  Is  CO/  6 


(56) 


Because  the  system  has  no  external  forces  applied,  it  remains  torque-free.  Thus,  Equation 
(28)  can  be  used  to  eliminate  ^c  co,  Q)i  term  and  arrive  at 


Ed  total  =  [Is  COs  + 1/  0)s') 


[It  total  ~  Is  total)  COs  —IsCJ 
It  total 


(57) 
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Noting  that  the  platform  nutational  frequency.  Equation  (33),  can  be  written  as 

V       1         r.'      [h(Os^Is   0)s\      ,    ,      {Js  total  -  h  to/a/)  6)5'  +  /j  O" 

A  =  Ao  -  G)y  = -(Os (58) 

V        h  total        I  h  total 

then  Equation  (57)  can  be  written  as 

Ed  total  =  [h  4  +  Is  ^s]  (-  A')  =  io  A  (-  X)  (59) 

Referring  to  Equation  (34),  the  nutation  angle  can  be  written  as 

"°^  ^  =  I ^i^ J  ^^^ 

Taking  the  derivative  and  comparing  it  to  the  rate  of  total  energy  dissipation.  Equation  (59), 
the  following  relationship  can  be  established 

_  0  sin  0  =  P-^^Hr^l  =  ^^^^  (61) 

Because  the  rotor  is  rigid,  all  energy  dissipation  will  occur  in  the  platform,  such  that 

—  =  |h|0sine  (62) 

A' 

Equations  (61)  and  (62)  can  be  rewritten  by  including  the  motor  torque,  Equation  (49),  as 

-|h|  d sin  d  =  Tr  +  //  0)/  =  -^  (63) 

A' 

Because  Tp/^  =  -  Trim  (action  -  reaction  pair),  the  final  relationship  is  written  as 

TpiM  =  —-^IsCOs'  (64) 

A' 

This  equation  describes  the  motor  torque  on  the  quasi-rigid  platform  as  a  function  of 
platform  variables  only.  It  can  be  seen  at  this  point  that  the  classical  analysis  can  be 
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achieved  by  assuming  no  torque  is  applied  by  the  motor,  resulting  in  TpiM  -  0  and 

^  =  -1/6)/  (65) 

This  analysis  can  now  be  performed  on  the  system  containing  a  quasi-rigid  rotor  and 
a  rigid  platform.  It  is  observed  that  the  selection  of  the  body  to  be  the  rotor  and  the  body  to 
be  the  platform  is  completely  arbitrary,  and  no  physical  distinction  exists  between  the  two 
bodies.  Therefore,  by  analogy,  the  equation  describing  the  torque  on  the  quasi-rigid  rotor 
as  a  function  of  rotor  variables  must  be 

T/f/M  =  ^  +  /.4  (66) 

A' 

Now  let  Tp/]^  and  Tp/j^  represent  the  motor  torques  on  the  platform  and  on  the  rotor 
respectively  for  the  system  containing  both  a  quasi-rigid  rotor  and  a  quasi-rigid  platform. 
Then  the  separation  axiom  would  require  the  following  two  conditions 

Tp/M  =  TpiM  =  — +  7/4'  (67) 

X' 

Trim  =  Trim  =  ^  +  />,  (68) 

A 

Once  again,  since  the  system  is  an  action  -  reaction  pair 


,* 


TpiM  +  TjiiM  =  0  (69) 


and  then 


^  +  ^  =  -[is^s^h'^s]  =  N^sin  e 


(70) 
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The  stability  critericMi  dictates  that  6  remain  zero  or  be  negative,  thus 

This  is  seen  to  be  the  existing  energy  sink  criterion.  Therefore,  despite  the  presence  of 
motor  torque,  the  existing  stability  criterion  still  applies.  Referring  to  Equations  (67)  and 
(68),  the  second  term  in  the  equations  would  represent  the  torque  applied  by  the  damping 
mechanism  to  the  platform  and  the  rotor  respectively.  Equations  (67)  and  (68)  state  that  the 
motor  torque  minus  the  damper  torque  will  equal  the  net  torque  of  the  platform  and  the 
rotor  respectively. 

In  determining  the  revised  energy-sink  stability  theory,  the  energy  dissipation 
equation  for  the  system  with  no  energy  contribution  fipom  the  motor.  Equation  (37),  must 
be  rewritten  to  account  for  the  motor  torque.  Thus 

E  =  Eototai  +  W  =  -/,Aa),-//A'(y/  (71) 

Substituting  in  Equations  (67)  and  (68),  one  arrives  at 


E  =  Eotoioi  +  W  =  Xl^-Tl 


RIM 
X 


+  ^'"^-^mi 


(72) 


-  Ed  +  Ed'  -  (A  Tffif^  +  X'  Tpn^) 
By  assigning  the  motor  torque  values  as 

Tm  =  T)iiM  =  -TpiM  (73) 

then  the  above  equation  can  be  rewritten  as 

E  =  £D/ora/  +  W^  =  Ed  +  ^d'  +  T-a/U'-A)  (74) 

Because  the  total  time  rate  of  change  of  energy  dissipation  equals  the  rate  of  change  of 
energy  dissipation  in  the  rotor  plus  that  in  the  platform,  the  rate  of  work  due  to  the  motor 
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torque  must  be 

W  =  TM(X'-X)  =  TMa  (75) 

Referring  to  the  example  from  Chapter  1,  the  prolate  dual-spin  communications 
satellite,  typically  the  rotor  nutation  frequency  would  be  negative  and  the  platform  nutation 
frequency  would  be  positive.  For  rotor  spin-up  Tm>0  and  the  motor  torque  would  add 
energy  to  the  system,  and  for  rotor  spin-down  Tm<0  and  the  motor  torque  would  remove 
energy  from  the  system.  For  the  arbitrary  system,  the  motor  torque,  T^,  and  the  sign  of 
the  term  X'  -X  (=  o),  would  determine  whether  the  motor  adds  or  removes  energy  from 
the  system. 

The  rates  of  change  of  the  energies  of  the  system  may  now  be  represented.  Rewriting 
Equation  (46)  to  determine  the  total  energy  dissipation  of  the  system 

EDtotai  =  E-W  (76) 

The  rate  of  work  due  to  the  motor  torque  can  be  expressed  by  combining  Equations  (75) 
and  (68)  (or  Equation  (67) )  to  arrive  at 

W  =  ^c+IsQ)s<y    =  -—  <T-// 0)/  a]  (77) 

A  \       A'  / 

The  rate  of  change  of  the  kinetic  energy  of  the  system  as  a  function  of  platform  core  energy 
and  system  parameters  was  determined  previously  as 

E  =  Ec+  Is  c(di/  +  ct)  +  Is  CO/  6  (52) 

Substituting  the  above  two  equations  into  Equation  (76),  the  expression  for  total  energy 
dissipation  of  the  system  may  now  be  written  as 

Ed  total  =  Ec+Is<^(o>/  +  CT)  +  IsCOs'cT-—a-IsQ)sO  (78) 

A 
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Because  O),'  +  6"  =  o),,  this  equation  can  be  further  simplified 

Ed  total  =  Ec'  +  IsC0/(T-^o  (79) 

A 

Comparing  this  equation,  representing  the  system  with  a  quasi-rigid  rotor  and  a  quasi-rigid 

platform,  with  Equation  (54),  representing  the  system  with  a  rigid  rotor  and  a  quasi-rigid 

platform,  reveals  an  additional  term,  — ^  o.  The  first  term  of  Equation  (79)  represents  the 

A 

rate  of  change  of  core  energy,  with  the  platform  as  the  core  body.  The  second  term  will 
account  for  the  change  in  energy  associated  with  a  change  in  the  relative  rotation  rate  of  the 
rotor  with  respect  to  the  platfomi.  The  final  tenn  accounts  for  the  energy  loss  due  to  the 
quasi-rigidity  of  the  rotor  that  is  not  represented  in  the  platform  core  body  expression. 

The  case  that  will  be  analyzed  is  that  which  occurs  when  there  is  a  constant  relative 
rotation  rate  between  the  platform  and  the  rotor.  For  the  remainder  of  the  analysis,  let 

<T  =  0  (80) 

and  the  total  rate  of  energy  dissipation  becomes 

Ed  total  =  Ec'-^^^  (81) 

A 

Further  simplification  can  be  achieved  by  noting  that  cr  =  X'-X  and  from  Equation  (37) 

that  Ed  total  =  Ed  +  Ed'.  Therefore 

Ed  +  Ed'  =  Ec  -  Ed  h— ^  (82) 

which  reduces  to 

^  +  ^  =  ^  (83) 

At         At  At 

A  similar  development  can  be  performed  using  rotor  core  energy  vice  platform  core  energy. 
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The  result  is 

^  +  ^  =  ^  (84) 

There  is  an  expected  symmetry  between  Equations  (83)  and  (84),  due  to  the  arbitrary 
assignment  of  one  body  of  the  system  as  the  rotor,  and  the  other  body  as  the  platform.  To 
confirm  the  results  of  Equations  (83)  and  (84),  one  must  prove  that  when  a  is  constant 


X        X' 


(85) 


(86) 


Taking  the  derivative  of  the  rotor  core  energy.  Equation  (48) 

^C  =  U  total  COt  6)t  +  Is  total  (Os  0)s 

=  ft  total  COt  6)t  +  Is  total  (o)/  +  c)  (w/  +  o) 

and  similarly,  for  the  platform  core  energy 

Ec     =  It  total  COt  0)t  +  Is  total  (Os   0)s'  (87) 

Substituting  Equation  (87)  into  Equation  (86)  and  noting  that  <T=  0, 

Ec  =  Ec'  + Is  total  <^o)/  (88) 

Multiplying  through  by  the  platform  nutation  firequency 

Ec  X'  =  Ec'  X  +  Is  total  X'  a  6)/  (89) 

and  recalling  that  X'  =  X+  o, 

EcX'  =  Ec'X  +Ec  a-^IstotaiX (ro>s'  (90) 

which  is  the  same  expression  as  Equation  (85)  if  it  can  be  proven  that 

Ec'^h  total  X'Ws'  =  0  (91) 

Eliminating  the  transverse  angular  velocity  of  Equation  (47)  by  substituting  in  Equation 
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) 


(26),  one  airives  at 

Ec    =  ^  J -^^  Is  total  CDs  (92) 

Taking  the  derivative 

,  ^       (/^  6)5  +  f/  (O/)  (is  Q>5  +  //  O)/)  ^  J  ...  /  .:.  /  //  total 


Ec 


and  simplifying 


If  total 


total  CO/ CO/ '-^^^  (93) 

*t  total 


t  /         \I?  COs  (Os  +  Is  1/  COs  d)/  +  Is  1/  0)/  C^,       ■/    0)/  6)sf  +  Is  total  It  total  CO/  6)/ 

Lc  = J — K^^) 

h  total 

Noting  that  Is  total  =  ^^  + 1/  and  recalling  that  O),  =  co/  because  cr  =  0,  the  above 
equation  can  be  reduced  to 

£.   /  ,  [is  (Os-^  1/  (O/  -  It  total  (O/  )  6)s  .OCA 

Ec    =  -Is  total J (95) 

't  total 

finally,  invoking  the  definition  of  the  platfonn  nutation  frequency,  one  arrives  at 

Ec    =  -  Is  total  ^'o)s  (96) 

Therefore,  Equation  (91)  holds  and  the  revised  stability  criterion  can  be  expressed  as 

(97) 

A  few  remarks  can  be  made  concerning  this  stability  criterion.  The  third  expression  is  the 
existing  energy  sink  stability  criterion.  This  criterion  must  equal  the  stability  criterion  as  a 
function  of  the  rotor  core  energy  which  must  equal  the  stability  criterion  as  a  function  of  the 
platform  core  energy.  It  can  be  seen  that  one  no  longer  needs  to  know  the  energy 
dissipation  rates  in  the  platform  and  in  the  rotor  to  determine  stability.  By  knowing  or 
postulating  the  core  energy  over  time,  the  stability  of  the  system  can  be  determined. 


Ec 

Ec' 

_Ed 

^Eo' 

< 

0 

X 

X 

A 

X 
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Continuing  with  the  prolate  dual-spinner  example  of  Chapter  1,  any  one  of  the  above 
expressions  apply.  For  a  stable  system,  the  rotor  core  energy  will  be  positive,  and  increase 
over  time.  Additionally,  the  rotor  nutation  frequency  will  be  negative,  resulting  in  a 
negative  expression  for  the  rotor  core  energy  stability  criterion.  The  platform  core  energy 
will  be  positive,  but  will  decrease  over  time.  The  platform  nutation  frequency  will  be 
ix)sitive,  resulting  in  a  negative  expression  for  the  platform  stability  criterion.  According  to 
Equation  (97),  the  platform  and  rotor  stability  criteria  must  equal  one  another.  From  the 
numerical  simulation  of  a  dual- spin  quasi-rigid  axisymmetric  system,  the  rotor  and  platform 
core  energies  as  a  function  of  time  will  be  determined  and  graphed. 

D.    NUTATIONAL  MOTION 

The  development  of  a  modified  expression  for  the  nutation  angle  as  a  function  of  time 
may  now  be  presented.  The  actual  nutation  angle  of  the  system  is  defined  as  6.  The 
nutation  angle  as  a  function  of  platform  core  energy  will  be  represented  by  rf,  and  the 
nutation  angle  as  a  function  of  rotor  core  energy  will  be  represented  by  7).  The  derivation  is 
similar  to  the  one  previously  given  in  this  chapter,  except  that  the  total  energy  of  the  system 
has  been  replaced  by  the  core  energy  of  the  system  to  eliminate  die  transverse  angular 
velocity,  fi),.  The  derivation  for  the  platform  core  energy  will  be  shown.  From  before,  the 
angular  momentum  of  the  dual-spin,  quasi-rigid,  axisymmetric  system  is 

h^   =   ^t^total  ^/^  +  (h  COs  +  /,'  (Osf   =   If^otal  CO,^  +  (ls(^+  h  total  (O/)  (26) 

Combining  this  with  the  platform  core  energy,  and  solving  for  the  platform  angular  velocity 
about  the  spin  axis 

h  total [h  total  -  It  total)  CO3'    -2  Is  total  h  O" (O^   +  //  (T^  +  2  Ec  U  total- hp-   =  0   (98) 


26 


and 


e>3'  =  7j—^, r±Vg(       '"^f       )-L- 

Vt  total  —  *s  total]  Vt  total  —  is  total)  's  total 


(99) 


where 


e'  =  2  E/  Is  totalis  -^i-^-h^{^f-^il  -^f^-^  (100) 

V    It  total'  ^It  totally        It  total'      \It  total' 

The  initial  conditions  at  r  =  0  will  determine  the  proper  sign  of  Equation  (99).  As  before, 
continuity  will  ensure  this  sign  for  all  f  >  0.  Substituting  Equation  (99)  into  Equation  (11) 

., (/.,»,^^3'+/.g|  ^  ^^,, ([,   ^""'f'    ]{[hc±  /e^l  (101) 

^  n  '  \[It  total  -  Is  total] "  / 


77    =  COS" 

The  time  rate  of  change  of  the  nutation  angle  will  detennine  the  stability.  Differentiating  the 
above  equation  results  in 


W  =  ~^^     1  /  h  total  \ 

Sinn  h  [±^j 


(102) 


A  stable  system  would  require  that  the  lower  sign  be  chosen  for  the  radical,  thus 

r,'  =  cos-^  l\.^uoml ^]i[/,cT-VG^)  (103) 

IL'/  total  ~  h  total!  "■  I 

When  this  sign  selection  is  applied  to  Equation  (52),  one  gets 

0>3'  (It  total  -  Is  total)  ^IsO^-^fQ'  (f^^l  (104) 

This  can  be  rewritten  as  the  well  established  dual-spin  stability  condition,  as  written  by  P. 
C.  Hughes  [Ref  6], 

(//  total  -  Is  total)  Oyi-^hO    <  0  (105) 

It  is  important  to  note  that  this  analysis  does  not  produce  a  contradiction  to  Equation  (41). 
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In  a  similar  manner,  the  modified  nutation  angle  can  be  derived  with  respect  to  the  rotor  as 

7,  =  cos-W[-— ^^-ifif ^IK^^'^-^]]  (^^) 

M't  total  ~  h  total]  "  / 


where 


Q  =  2EcIs 


roro/ (l  -  f^^l  - /I^  f^^-^l  (l  -  t^-^l  +  (l^ 
»    *t  total'  y't  total' ^        't  total'      v// 


total' 


(107) 


Therefore,  if  one  is  given  the  core  energy  or  the  postulated  core  energy  as  a  function  of 
time,  the  nutational  motion  can  then  be  determined.  This  leads  to  an  extremely  important 
conclusion.  By  determining  a  sufficientiy  accurate  postulation  of  the  core  energy  of  a  dual- 
spin,  quasi-rigid,  axisymmetric  spacecraft  over  time,  one  can  predict  the  nutational 
behavior  and  the  stability  of  that  spacecraft.  Numerical  simulation  of  such  a  system  will 
confinn  this  conclusion. 
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in.  DYNAMICAL  EQUATIONS 

A.    MINGORI  DUAL-SPIN  SYSTEM 

A  specific  mcxiel  is  required  to  validate  the  mcxlified  energy-sink  stability  theory.  The 
development  of  the  previous  chapters  was  completely  general  in  nature  and  applies  to  any 
dual-spin  spacecraft  with  an  arbitrary  damping  mechanism.  D.  L.  Mingori's  dual-spin 
system  [Ref.  3]  provided  the  needed  model  required  to  validate  the  proposed  stability 
theory.  Figure  (3).  Additionally,  the  complete  non-linear  equations  of  motion  were 
presented  by  Mingori;  however,  the  expressions  for  the  important  parameters  in  attitude 
dynamics,  angular  momentum  and  kinetic  and  total  energy,  were  not  presented  in  his  paper 
and  had  to  be  derived  before  any  analysis  could  be  performed. 

The  Mingori  system  is  comprised  of  two  symmetric  rigid  bodies,  the  lower  which 
shall  be  defined  as  the  rotor,  and  the  upper  body  shall  be  defined  as  the  platform.  By 
convention,  all  terms  referring  to  the  platform  will  be  the  same  notation  as  that  of  the  rotor, 
except  that  they  will  carry  the  prime  mark.  Both  the  rotor  and  the  platform  have  coordinate 
axes  fixed  to  the  body  and  located  at  the  respective  centers  of  mass.  The  distance  between 
the  centers  of  mass  is  specified  by  L.  The  entire  spacecraft  has  coordinate  axes  fixed  to 
the  spacecraft  center  of  mass,  denoted  by  the  double  prime  coordinate  axes,  and  rotating  in 
the  same  manner  as  the  rotor  coordinate  axes.  The  coordinate  axes  bs,  bs',  and  b3"arc  all 
coUinear.  The  spacecraft  center  of  mass  will  vary  along  the  bs  axis  as  the  point  masses  in 
the  rotor  and  the  platform  oscillate.  The  only  relative  motion  of  the  platform  with  respect  to 
the  rotor  is  angular  rotation  about  the  bs  axis.  A  motor  driven  by  a  control  system 
maintains  a  constant  relative  rotation  rate  o.  The  angle  between  a  line  parallel  to  bi  and 
passing  through  the  platform  center  of  mass,  and  bi'  is  represented  by  v^  =  a/ . 
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»»3.b'3.b"3 


Mingori  Dual-Spin  Quasi-Rigid  Axisymmetric  System 
Figure  (3) 
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Each  body  contains  a  mass-spring-dashpot  mechanism.  An  actual  spacecraft 
undergoes  damping  from  various  mechanisms.  Undesired  damping  can  occur  due  to 
structural  elasticity  and  liquid  propellant  slosh.  Nutation  dampers  are  incorporated  into 
spacecraft  to  improve  the  nutational  motion.  The  dual-spin  axisymmetric  system  with 
mass-spring-dashpot  dampers  cannot  accurately  model  an  actual  spacecraft,  but  is  used  to 
illustrate  and  validate  the  theory.  A  description  of  the  rotor  is  as  follows.  The  mass- 
spring-dashpot  mechanism  can  be  nxxlelled  by  a  particle  mass  mi  attached  to  a  spring  with 
constant  k  inside  a  tube  filled  with  viscous  fluid  with  damping  coefficient  c.  The  motion 
of  the  particle  is  constrained  parallel  to  the  bs  axis  only.  At  rest,  the  particle  mass  lies 
along  the  rotor's  bi  coordinate  axis,  at  a  distance  a  from  the  rotor's  center  of  mass.  Three 
balancing  masses,  mi,  m^,  and  m4,  each  equal  to  the  mass  of  the  mass-spring-dashpot 
mechanism,  are  rigidly  fixed  a  distance  a  on  the  b2,  -bi,  and  -b2  axes.  These  masses 
maintain  the  axisymmetry  of  the  system  about  the  bs  axis.  Displacement  of  the  particle 
mass  mi  is  denoted  by  the  variable  z.  A  simplification  in  this  paper  of  the  Mingori  dual- 
spinner  system  is  the  assumption  of  a  massless  spring-dashpot  system.  Thus  the  particle 
of  the  mass-spring-dashpot  system,  mi,  is  the  same  mass  as  the  corresponding  three  other 
balancing  masses,  mi,  ms,  and  m^.  The  platform  can  be  described  in  a  similar  manner, 
with  all  notation  nwdified  with  the  prime  mark. 

The  dual-spin  system  center  of  mass  coordinate  axes  rotates  at  the  same  rate  as  the 
rotor  coordinate  axes,  but  is  located  along  the  bs  axis  at  a  distance  Zcm  ^^  the  rotor  center 
of  mass.  This  distance  will  vary  as  the  particle  masses  are  displaced.  Relating  the  dual- 
spin  system's  coordinate  axes  to  the  rotor  coordinate  axes  was  arbitrary.  Equivalent  results 
would  be  achieved  by  selecting  the  platform  coordinate  axes  instead. 
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B.    ANGULAR  MOMENTUM 

The  angular  momentum  of  the  Mingori  dual-spin,  quasi-rigid,  axisymmetric  system  is 
now  derived  from  first  principles. 

The  angular  momentum  (moment  of  momentum)  of  a  particle  of  mass  m,-  located  in 
body  B  is  defined  as 

hi  =  riXmiV  (108) 

where  h,  is  the  angular  momentum  of  the  /  th  particle  with  respect  to  the  system  center  of 
mass,  Fj  specifies  the  location  of  the  /  th  particle  with  respect  to  the  system  center  of  mass, 
andV  =  4^ 


.  is  the  absolute  velocity  of  the  /  th  particle  with  respect  to  the  inertial 

at    N 

reference  frame  N.  Expressing  the  absolute  velocity  in  the  system  reference  frame 

V  =    V     +  Fi  +    CO   X  r,-  (109) 

where  B  represents  the  reference  coordinate  axes  of  the  system.  For  the  following 
derivation,  all  displacements  and  velocities  are  referenced  with  respect  to  the  rotor  (b,) 
coordinate  axes.  The  angular  momentum  vector  is  rewritten  as 

h,  =  r,  X  m,  \  V     +  r,  +    co   x  r,j  (1 10) 

Applying  this  equation  to  particle  1  with  mass  mi  and  position 

ri  =  a  bi  +  0  b2  +  (z-z^mjbs  (111) 
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one  amves  at 


(112) 


hi  =  mi 


(Z'Zcrrf  0  -a(z-Zcm) 

0  a^^(z-zcmf  0 


.  -a(z-Zcm) 


0 


fi)l 

©3 


0  -(2-Zcm)  0 

(Z-Zcm)  0  -a 

\L      0  fl  0 


Vany      + 


0 


-fl(z-Zcm) 

0  J/ 


'11 

ibsJ 


For  the  rigid  body.  Equation  (1 10)  is  applied  to  a  differential  volume  at  location 

r^Af  =  "  bi  +  V  b2  +  (w  -  Zcn^  b3  (1 13) 


Integrating,  one  arrives  at 


/ 


hA/  = 


h  +M  Z}m  0 

0  h+Mz}, 

0  0 


cm 


0 
0 


Ct)i 
0)2 

a>3 


0 

-A^Zcm 

\L    0 


MZcm 

0 
0 


0 
0 
0 


cm  X 


cm  y 

cm  z 


"1 

b2 
b3 


(114) 
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Similar  calculations  are  performed  for  the  remaining  seven  point  masses  and  the  platform. 
For  the  rotor,  one  arrives  at 

2cm 
2em 


m 


I    Z^-2zZcm    / 


-maz 


-maz 


h  +  M  z}m  + 

'"\  z^-2zzcm  I 
0 

0 

-maz 

0 

(Af  +  4m)(2c«) 
-  m  z 


0 

{M  +  4m){zcm) 
-  mz 


h  +  4  nup- 


0 
0 
0 


0)1 


a>2 


0)3 


\ 


/       \ 

Vo», 

V'o-y 

»' emx 

b2 
lb3, 
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Similarly,  fOT  the  platform 


hM' 


/2'+Af'(/-2c-.f  + 

'^  U''  +  2z'(/-Zc.) 
-  m'  a'  z'  cos  (<t  /)  -  m'  a'  z'  sin  [a  t) 


-  m'  a'  z'  cos  (a  f) 


-  m'  a'  z'  sin  (<y  /) 


h'  +  Am'a'^ 


(0\ 


'    a>2 


fl>3 


-  m'  fl'  <Tcos(a  /)  z'  +  m'  a'  sin(a  /)  z' 

-m  a'  a  sin(cT  /)  z'  -  m'  a'  cos(<r  i)  z 

(/3'  +  4m'fl'^)cr 


-(A/'  +  4m'){/-Zc«) 


(A/'  +  4m')(/-Zc«) 


(116) 


The  angular  momentum  equations  have  terms  corresponding  to  the  velocity  of  the 
center  of  mass.  Angular  momentum,  when  taken  about  the  center  of  mass  of  a  system,  by 
definition  must  be  independent  of  the  translation  of  the  system's  center  of  mass  (with 
respect  to  an  inertial  reference  frame).  To  verify  that  this  occurs  in  the  above  equations,  the 
equation  for  the  position  of  the  center  of  mass  is  substituted  into  the  center  of  mass  velocity 
terms,  and  then  these  equations  are  equated.  If  they  have  the  same  magnitude  but  opposite 
sign,  they  will  cancel  each  other  out  and  it  will  be  proven  that  the  system  angular 
momentum  is  independent  of  the  translation  of  the  system  center  of  mass. 
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One  must  prove 

Looking  at  the  bi  components,  one  finds 

-[(A/'  +  4m')(-/  +  z,„)-m'z']=[(Af  +  4m)z,«-m2]  (118) 

Rewriting  this,  one  obtains 

(M'  +  4m')l  +mz  +  mz'  =  (M  +  4m  +  M'-i-4m')zcm  (119) 

The  center  of  mass  with  reference  to  the  b  i ,  b2,  b3  coordinate  axes  is 

y,  Mi  Zi  +  m,-  Zi  /  X 

i  mz+  A/'  +  4m'  /  -^  m' z' 

2cm  =  -^^= =  ^^ (120) 

2^Mi  +  nti  M  +  4m  +A/'  +  4  m' 

< 

Substituting  Equation  (120)  into  Equation  (1 19)  results  in 


(.J,     A     At                  /   /  ?  ^^     -.          MMf     A     Aimz  +  yM' +  4m')l  +m' z' 
[M  +4m)l  +mz  +  m  z    =  (A/  +  4m +  A/  +  4m  )  ^^ 

\      M  +  4m  +M'  +  4m'      ^ 


(121) 


This  reduces  to 

(M'  +  4m)l  +mz  +  mz'  =  (M'  +  4m)l  +mz  +  mz'  (122) 

Therefore,  the  angular  momentum  along  the  bi  axis  is  independent  of  the  velocity  of  the 
center  of  mass.  Verifying  this  along  the  other  axes  can  be  done  in  a  similar  manner.  This 
proves  that  the  angular  momentum  of  the  entire  system  about  the  system  center  of  mass  is 
independent  of  the  translation  of  the  system  center  of  mass.  The  angular  momentum  of  the 
entire  system  can  be  written  as 
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where 


hM*Am  = 


1 

/i+Afz2,+ 
^(2a2  +  4ziL+\ 
'"I  z^-2zzc.  1 

0 

-maz 

0 

/2  +  A/22.+ 

'"1  22-2rzc-.  / 

0 

-maz 

0 

/3  +  4  ma^ 

0 

1 

+ 

-maz 

\ 

0      . 

and 


/r+Af'(/-zc^f+ 

^,|2fl''  +  4(/-Zc„f  + 
z'2+2z'(/-2cm) 


0 


/2'+M'(/-Zc«f  + 

^,/2fl'^  +  4(/-Zc».f  + 
z'2+2z'(/-2.^) 


6>1 


0)2 


0)3 


-  m'  a'  I  cos  (<7  r)  -  m'  a'  z'  sin  (<t  f) 


-  m'  a'  z'  cos  (<t  /) 


-  m'  a'  z'  sin  (<t  /) 


/3'  +  4m'fl'^ 


-m  a'  o  cos(a  r)  z'  +  m'  a'  sin{<T  r)  z' 
-m  a'  a  sin(cT  /)  z'  -  m'  a'  cos(a  /)  z' 

(/3'  +  4m'a'V 


bi 

b2 

ib3 


0)1 


0)2 


0)3 


(124) 


(125) 


b2  / 
ib3| 
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C.    ENERGY 

The  kinetic  energy,  E,  and  the  total  energy,  Etotai^  are  now  derived.    From  first 
principles  the  kinetic  energy  of  a  differential  particle  is  written  as 

^^  =  ^hv'"pdm  =  i-((^v'").  (^V'^))^  (126) 

where  the  absolute  velocity  is  defined  from  before 


V     =    V     +rrfm+    CO   xrrf;„ 


(109) 


The  kinetic  energy  of  the  Mingori  system  is  detennined  by  integrating  the  differential 
kinetic  energy  over  the  rotor  and  over  the  platform,  and  summing  the  differential  kinetic 
energy  equation  over  the  eight  point  masses  to  arrive  at 

f  ((V")(v)idm+ii{(v)(v))'".* 

A—  1=1'' 

(127) 

;n^*^  r  =  1 


E  =  ^ 
^        2 


Performing  the  steps  on  particle  /n  i ,  the  following  is  obtained 


(128) 


^^^2"\.f^J2(V'".n)^2(V'".(V'"xn))+2(h-(^5^'"xn)) 


Substituting  in  the  values,  one  arrives  at 


2  .   \ 


El  =^mi 


'y}mx+  Vcm  >  +  V2„  2  +  (z  -  Zcmf  +  {z  -  Zcmf  CO^  + 

fl2  G)J  +  (z  -  Zcmf  a)\-2a{z-  Zcm)  G)l  0>3  + 

fl2G)|  +  2(z  Vcmz-icm  Vcmz)  + 

I    Z  fi>l  Vcm  >  +  2cm  (0\Vcmy-aa)2  Vcm  z   ] 
2(-a'zQ>2  +  aZcmO>2) 


(129) 
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This  operation  is  then  performed  over  the  other  three  panicles  in  the  rotor.  Performing 
these  same  steps  on  the  body  of  the  rotor, 

(130) 


Em  =  ^ 


(V"  •  V)  +  (h  •  ri)+  ((^S^'^x  ri  )  .  (^5-x  ri) )  + 


2|^^l2(V"-rJ+2(V'".(V"xn))^2(n.(V"xrO) 


dM 


Substituting  in  values,  and  specifying  the  position  vector  of  the  differential  particle  as 

r^A/  =  "  bi  +  V  b2  +  (w  -  Zcm)  bs  (131) 

then 

(132) 


I 


i„-i 


I 


Z^m  6J2  +  v2  fi)|  -  2  W  V    0)2  G)3  +  2  V  Zcm  Q>2  fi>3  +  ^2  6)?  + 
z2„  Ojf  +  W2  G)|  -  2  W   Zcm  ©?  -  2  U  W  0)1  0)3  +  2  M  Zcm  CO]  0)3  +   1 

^  v2  0)?-2vM0)i0)2  +  M2a,|  +  2(-ZcmVcmz)+  r^ 

^Iw  0)2  Vcmx-  2cm  0>2  Vcmx"  V  0)3  V^mx"  ^  0)i  V^m  y  +] 

2  (-  V  Zcm  0)1  +  M  Zcm   6)2) 


Integrating  over  the  limits  of  the  symmetrical  body,  the  kinetic  energy  of  the  rotor  becomes 

,2 


Em  =  ^M 


IVL  X  +  V?m  >  +  V?m  z  +  z|m  +  Z?m(a>?  +  fl>|)  A 

h  fi>?  +  /2  fi>2  +  ^3  ^3  - 
\  2  Zcm  fi>2  Vcmx  +  2  Zcm  fi>l  ^cm  j,  -  2  Zcm  Vcm  zl 


(133) 


These  same  steps  are  performed  on  the  platform  and  on  the  remaining  seven  particles. 
These  equations  are  then  combined  and  simpliHed. 
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The  total  kinetic  energy  for  the  Mingori  system  is 

(134) 

^(z^  +z^\C0i+e02]-2az  CDi  0)3)  + 

^ \z    +z'  \a>\  +  0)2) -la' z' cos(ot) a)i  0)3  +  2  j' z' sin(cTr) 0)2  g>3)  + 

m(-a'z  (02)  + 
£  =   I  m'(-a'2'cos(crr)G)2  +  fl'z'sin((Tf)fi)i-fl'2'<Tcos(crr)a)i-a'z'<Tsin(crr)G)2) 

1  (fi)f  +  0)2 ) (mz  +  m'  z'f  -(mz  +  m'  zf 


2  Mt 

..  (a/'  +  4  m'f 

+ 


lU^a)|){A/-H4m)/2'^'-'^'"'^' 


1  (o)?  +  fi)|)  (m'  +  4  m')  /2  (^If^^f  - 

(    9        9\(m  z  +  m' z' +  (A/' +  4mj/| 
m  z  [(o\  +  0)2')  ( -j^^ ) 

/   ,i    1        o\\mz-¥m' z' -{M -^Am)! 
mz  (0)1+ 0)2) 


Mj 


The  total  energy  of  the  Mingori  system  can  be  determined  by  adding  the  kinetic 
energy  of  the  system  to  the  potential  energy  of  the  system.  The  only  potential  energy  of  the 
system  is  the  energy  stored  in  the  springs  of  the  mass-spring-dashpot  damper.  From 
Hooke's  law,  the  system  potential  energy  is  written  as 

C/(z)  =  i/:z2  +  iit'z'^  (135) 

The  total  energy  of  the  Mingori  system  is  then 

Etotai  ^  E+U  =  E  +  ^kz^  +  ^k'z'^  (136) 
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D.    MINGORPS  EQUATIONS  OF  MOTION 

Mingori's  equations  of  motion  [Ref  3]  arc  written  as  follows 

A  Q)i-(A-C)o)2  cos-Js' CC02  +  2  Afr  f  f  fi)i  +  A/r  f  (6)1  +  0)2  0)3)  + 

-2  (f  +  l2)[z  (6)1  -  6)2  G>3)  +  i  CDi]  + 


m  \z 


-2  CC0i+2zQ)i  +Z  (6)1-0)2  G)3)-fl(d)3  +  £iDi6>2)j  ) 

-  2  (C- /i)[z' (d)i  -  6)26)3)  +  i' fi>i]  + 

-2  f  6)1  +  2z'6>i  +  2'(6)i-6)2^)-a'c»SV'^(6)3  +  6)16)2)]+  /  =  0 

i/smY  (z'  +  z' [(^  +  <^ - ^2 J I 

A  6)2  -  (C  -  A)  6)1  6)3  -  /s'  <T  6)1  +  2  Afr  f  f  6)2  +  A/j-  f  (6)2  +  6)1  6)3)  + 
I  -2  (f+ 72)  [2  (6)2 +  6)1  6)3) +  26)2]+  I 

|z[-2  ^6)2+226)2  +  2(6)2  +  6)1  6)3)J-fl[2  +  2(6)3 -6)f)]| 
-  2  (C-  /l)[z'  (6)2  +  6)16)3)  +  2'  6)2]  + 


m 


2' [- 2  f  6)2 +  2  2' 6)2 +  2' (6)2 +  6)16)3)]-  / 

fl'cosv^{2'  +  2' [(6)3  +  cf-0)i] )- a' sin v^ 2' (6)3 -6)1  6)2)/ 

C  6)3  -  m  fl  [2  2  6)1  +  2  {6)1  -  6)2  6)3)  ]  - 
m'  a'  sin  \fr[2  z  6)2  +  2'  (6)2  +  6)1  6)3)  J  - 
m  a'  cos  v^[2  z  (0\  +  z  (q)\  -  6)2  6)3)  J  =  0 

m(l  -p)z-m  pz  -ma(6)2-(0\  6)3)- 

m  (6)f  +  6)2 )  [2  ( 1  -  p)  -  /2  -  p'  2']  +  c  2  +  /:  2  =  0 

-mp'2  +  m'(l -p')2'  + 

m'a'{sinv/'[6)i  +  6)2(6)3 +  2  <t)]- cos ^[6)2- 6)1  (6)3  +  20)] ) 

m'(6)f  +  6)|)[2'(l-p')  +  /i-p2]  +  c'2'  +  i(:'z'  =  0 


=  0 


(137) 
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In  the  above  equations,  the  following  relationships  are  used 


V  s 


p  =  JO- 


Mt 


^      Mt 


(138) 


C  =  (pz  +  p'/)  Y^at 

A  =  /i+/i'  +  2mfcfl2  +  2mfc'a'^  +  {Af'  +  4mfc')(l-v)/2 

2 

C  =  /3  +  /s'  +  4  mfr  a2  +  4  ^^'  ^' 


73'  =  /3'  +  4mfe'fl'^ 

Chapter  IV  discusses  how  these  equations  are  adapted  for  use  in  the  numerical  integration 
routine. 


E.    NUTATION  ANGLE 

The  nutation  angle  of  the  dual-spin  system  can  be  detennined  by  substituting 
Equations  (124)  and  (125)  into  Equation  (34)  to  arrive  at 

(139) 


0  =  cos-M 


IN 


M  =  ..e-i^!lli!L  = 


cos 


|h| 


cos 


-1 


^  (-m  fl  z-m' a' z' cos(  (J  r))  0)1  +' 

(-m'  a'  z'  sin{at))o>2  + 

(/3  +  4  m  fl2  +  /j'  +  4  m'  fl'  f  0)3 

+  (/3'  +  4m'a'% 


\ 


|h| 
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F.    CORE  ENERGY 

The  platform  core  energy  and  the  rotor  core  energy  were  defined  in  Equations  (47) 
and  (48).  The  nutation  angle  as  a  function  of  core  energy  was  then  developed  and  is  stated 
in  Equations  (101)  and  (106).  These  equation  can  be  used  to  predict  the  nutation  angle  of 
the  system  as  a  function  of  time.  Because  the  core  energy  as  a  function  of  time  is  not 
available  for  the  prediction,  a  postulated  core  energy  must  be  developed.  The  initial  value 
of  the  actual  core  energy  and  the  postulated  core  energy  must  match  and  is  dictated  by  the 
initial  conditions.  Parameters  of  the  postulated  core  energy  must  be  selected  to  accurately 
model  the  actual  core  energy,  to  include  the  final  energy  state  and  the  rate  at  which  the  it 
approaches  the  final  energy  state.  Two  models  were  considered. 

1.  Exponential  Core  Energy  Model 

An  exponential  representation  of  rotor  and  platform  core  energy  as  a  function  of 
time  are  expressed  as  follows 

Ec  postulated  {t)  =  (Ec^  -  Ec  final)  C^" '"')  +  £c  final  ( 1 40) 

Ec  postulated  (t)  =  [ECo  "  Ec  final)  ^^~''''  +  Ec  final  (141 ) 

The  initial  core  energies,  Ecc  Eqq,  are  determined  by  the  initial  conditions  of  the  system. 
The  final  core  energies,  Ec  fmah  Ec  final  y  as  well  as  the  exponential  factors  r,  /,  must  be 
selected  The  methodology  for  selecting  these  values  is  explained  in  Chapter  V,  Computer 
Analysis. 

2.  Verhulst  Logistic  Core  Energy  Model 

The  exponential  model,  as  will  be  explained  in  Chapter  V,  has  excellent 
agreement  for  stable  conditions,  but  performs  poorly  for  the  unstable  conditions.  As  an 
altemative  model,  the  following  first  order  differential  model  is  selected  for  the  rotor  and 
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platform  rcsi>ectively. 


^C  postuUaedi^)  - 


{Eco){Ec  final) 


Ec  postulated  {t)  = 


ECo +  {EC  final -Eco)^^-''^ 
\ECo)(Ec final  } 


(142) 


(143) 


ECo'AEc  final' -Ec:)ci-''^) 

This  type  of  equation  was  first  introduced  by  P.  F.  Vcrfiulst  to  model  human  and 
other  populations  [Ref  7].  It  is  often  referred  to  as  the  Verhulst  equation  or  the  logistic 
equation.  Although  population  dynamics  appears  to  be  unrelated  to  the  stability  of  a  dual- 
spin  system,  the  behavior  over  time  of  this  equation  compared  to  the  stable  and  unstable 
systems  provides  some  insight.  Figure  (4)  shows  the  Verhulst  equation  versus  time  for 
varying  initial  conditions. 
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Figure  (4)   Verhulst  Logistic  Equation  Versus  Time 
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It  can  be  noted  that  if  the  initial  value  of  the  variable  (in  our  case  core  energy)  is 
greater  than  the  equilibrium  value,  it  will  approach  the  equilibrium  value  in  a  manner  similar 
to  that  of  an  exponential  decay.  For  initial  values  that  are  less  than  but  within  one  half  of 
the  equilibrium  value,  the  variable  will  asymptotically  approach  the  equilibrium  value.  If 
the  initial  value  is  less  than  one  half  of  the  initial  value,  the  variable  over  time  has  a 
somewhat  different  shape  as  it  approaches  equilibrium.  Initially  the  slope  is  very  small,  but 
increases  to  a  maximum  at  about  the  one-half  of  the  equilibrium  position.  After  this 
inflection  point  the  variable  approaches  equilibrium  asyn:q)totically.  The  value  of  ^o  equals 
AT  is  an  asymptotically  stable  equilibrium.  A  value  of  ^o  equals  zero  is  an  unstable 
equilibrium.  If  the  value  of  No  is  slightly  greater  than  zero,  then  N(t)  will,  as  t  -*«», 
achieve  the  stable  equilibrium  of  K.  These  types  of  curves  will  be  utilized  to  describe  both 
the  stable  and  unstable  dual-spin  system.  The  initial  core  energies,  £co»  Eco,  would  be 
determined  by  the  initial  conditions  of  the  system.  The  final  core  energies,  Ec  nnah 
Ec  final,  niust  be  known  or  a  best  estimate  used.  The  exponential  factors,  r,  /,  must  be 
determined  experimentally,  and  are  a  function  of  the  system  parameters  and  initial 
conditions.  Chapter  V  provides  additional  details  on  the  selection  process. 
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IV.  NUMERICAL  SIMULATION 


A.    NUMERICAL  SIMULATION  EQUATIONS  OF  MOTION 

The  equations  for  the  dynamical  quantities  needed  for  analysis  of  the  dual-spin 
system  were  previously  developed.  Manipulation  was  required  to  make  these  equations 
suitable  for  numerical  analysis,  and  is  described  below.  The  computer  code  is  included  as 
Appendix  B. 

1.  Mingori's  Equations  of  Motion 

The  equations  of  motion  for  the  dual-spin  system  are  listed  as  Equation  set  (137) 
and  (138).  The  Runge-Kutta  numerical  integration  routine  necessitated  that  the  equations 
of  motion  be  a  set  of  first  order  differential  equations.  Mathematical  manipulation  was 
required  to  put  them  in  this  form.  Variables  representing  groups  of  terms  (i4„  Bi,  Ci,  etc) 
were  introduced  to  simplify  the  manipulations  of  the  equations  and  provide  a  suitable 
format  for  incorporating  the  equations  into  the  computer  program.  Mingori's  equations  can 
be  expressed  in  terms  of  these  variables  and  the  five  time-dependent  variables  of  motion 

A26  fih  +  ^27  fi>3  +  Ma  2'  +  i428  =  0 
B25  Q>2  +  B23  6>i  +  Bi3Z  +  B2\  z'  +  ^26  =  0 

Cioa)i+C5G)2  +  CG>3  +  Cii=0  (144) 

Di  Z  +  £>2  z'  +  I>3  G)2  +  Z>8  =  0 

£1  2  +  £2  z'  +  £3  Oh  +  ^5  6)2  +  f  10  =  0 

where  the  variables  /I,,  Bi,  Cj,  Z)„  £,,  Fj,  and  Z,  are  defined  in  the  Notation  section. 
Clearly,  these  equations  are  highly  coupled.  Through  a  series  of  manipulations,  including 
substitution  and  combining  sets  of  equations  to  eliminate  common  variables,  one  can  arrive 
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at  a  series  of  first  order  differential  equations  suitable  for  numerical  integration  techniques 


^=^ 

dt 

dL  -  2  - 

^-e  -^*g 

J   —  *■  — 
dt 

--^: ,  ^^<: 

dz'  _ 

-3^g  ^-g 

dt    ' 

i^.^ 

(145) 


^-g     ^-g 

(C5  ^13\  5  ^  fCio  i424  ^  CVB2l\  5'  ^.  ^10^28  ^  C5  B26  _  Cu_ 

dO)2_  ^  ^  ^  \CB25I       \CA26        C  B25  f         CA26C  B25       C 
dt  I     Cio  A21     C5  ^23 

C>426         Cfi25 

"f  >^26  >126  ^26 

d(Ol    =   ^  =       ^23  ^      ^13  =      ^21  =^      ^26 
^  B25  ^25         ^25  ^25 

The  equations  remain  coupled,  so  the  sequence  in  which  the  equations  are  numerically 
integrated  is  important.  Because  the  equations  for  z  and  the  equation  for  z'  are  functions  of 
z  and  z\  z  and  z'  must  be  integrated  first.  Similarly,  z  and  z'  must  be  integrated  before  ci>3, 

G)3  before  g>\,  and  6>\  before  Wi.  In  the  computer  program,  the  seven  variables  that 
describe  the  motion  are  stored  in  a  seven  column  matrix  where  co\  =  yll][j],  a>2  =  y[2]  [jl, 

fi)3  =  y[31  Ij],  z  =  y[4]  [j],  z  =  yI51  [jl,  z'  =  yI6]  Dl,  and  z'  =  y[7]  Dl.  The  index  j  identifies  the 
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matrix  location  for  each  set  of  saved  variables  over  time.  All  other  needed  quantities  can  be 
calculated  by  using  the  seven  time  dependent  variables  describing  the  motion,  and  the 
parameters  of  the  dual-spin  system.  These  additional  dynamical  quantities  are  stored  as  one 
dimensional  vectors  in  the  computer. 

2.  Angular  Momentum 

Angular  momentum  was  derived  in  Equations  (123),  (124),  and  (125).  The 
angular  momentum  is  a  vector,  but  for  verification  of  the  conservation  of  angular 
momentum,  only  the  magnitude  is  required.  Therefore,  the  three  vectorial  components  for 
both  the  rotor  and  the  platform  are  computed  individually  as  HI,  H2,  H3,  H1p,  H2p,  H3p, 
and  then  the  magnitude  of  the  angular  momentum,  h[j],  is  determined  and  stored. 

3.  Nutation  Angle 

The  nutation  angle  is  determined  using  the  previously  determined  relation 

Since  it  is  a  function  of  angular  momentum,  it  can  now  be  calculated  and  stored  as 
theta  [j]. 

4.  Energy 

The  total  energy  of  the  system.  Equation  (136),  is  determined  by  the  sum  of  the 
kinetic  energy,  as  derived  in  Equation  (134),  and  the  potential  energy  of  the  particle  masses 
of  the  mass-spring-dashpot  damper  system,  Equation  (135).  It  is  written  as 

Etotai  =  E-^U  =  E  +  ^kz^  +  ^k'z'^  (136) 

In  the  computer  program  the  kinetic  energy,  ke  [jl,  is  the  sum  of  the  constituents  of 
Equation  (134),  identified  as  T1  through  T13.  The  potential  energy  is  not  explicitiy  stated 
in  the  program,  but  instead  is  included  in  the  calculation  of  the  total  energy,  etotai  [j]. 
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5.  Core  Energy  and  Postulated  Core  Energy 

The  platfonn  and  rotor  core  energy  of  the  dual-spin  system  is  defined  in 
Equations  (47)  and  (48)  respectively.  All  of  the  variables  are  available,  such  that  in  the 
computer  program  the  energies  are  readily  calculated  as  Ecp  [fl  and  Ec  \j].  The  postulated 
nutation  angle  may  now  be  computed.  Equation  (101)  indicates  a  sign  must  be  selected 
depending  on  whether  the  system  is  stable  or  unstable.  The  calculation  of  nutation  angle 
over  time  in  Equation  (34)  will  indicate  stability.  The  computer  program  incorporates 
conditional  statements  to  assign  the  comet  sign  in  the  postulated  nutation  angle  expression. 
Equations  (101)  and  (106).  In  the  derivation  leading  up  to  the  postulated  nutation  angle,  o 
has  been  defined  as  the  relative  rotation  rate  of  the  rotor  with  respect  to  the  platform,  and  is 
a  positive  value.  This  provides  a  positive  contribution  to  the  angular  momentum  vector, 
thus  providing  stability  to  the  system.  In  Mingori's  equations  of  motion,  the  reference 
coordinate  axes  are  fixed  to  the  rotor,  resulting  in  the  relative  rotation  of  the  platform  in  the 
counter-clockwise,  or  negative,  direction.  The  postulated  nutation  angle  equations  were 
developed  using  the  former  reference  frame.  To  compensate  for  the  difference  in  the 
reference  frames,  the  postulated  nutation  angle  equations  in  the  computer  program  have  o 
replaced  by  -a. 

The  postulated  core  energy  as  a  function  of  time  is  computed  for  the  platform  and 
the  rotor,  using  Equations  (140)  and  (141)  for  the  exponential  model  and  Equations  (142) 
and  (143)  for  the  Verhulst  model.  A  postulated  core  energy  model  that  accurately  models 
the  actual  core  energy  required  several  iterations  to  find  the  proper  values  of  the  exponential 
factor. 

Once  the  postulated  core  energy  models  are  computed,  the  Q'  and  Q  values,  as 
defined  in  Equations  (100)  and  (107),  are  determined.  The  postulated  nutation  angles. 
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etah  [j],  and  eta  [j],  are  computed  using  Equations  (101)  and  (106)  with  the  proper  sign 
previously  detennined. 

The  revised  stability  criterion,  Equation  (97),  requires  the  time  rate  of  change  of 
the  core  energy  and  the  nutation  frequency.  Since  both  of  these  parameters  may  vary  with 
time,  an  estimate  is  made  to  determine  if  the  modified  stability  criteria  correctly  predicts 
stability  or  instability.  The  time  rate  of  change  of  the  core  energy  is  approximated  by  taking 
the  difference  of  the  final  and  the  initial  core  energy  values,  and  dividing  by  the  time  of  the 
simulation.  This  quantity  is  then  divided  by  the  initial  nutation  frequency.  This  stability 
quantity  is  computed  for  both  the  rotor  and  the  core,  but  is  indicated  in  Table  (1),  Summary 
of  Analyzed  Cases,  simply  as  a  negative  or  positive  value. 

B.    COMPUTER  PROGRAM 

Essential  to  the  verification  of  the  stability  criteria  is  the  computer  program  that 
integrates  the  dual-spin  system  equations  of  motion,  calculates  other  dynamical  quantities 
of  motion,  and  graphs  the  results.  The  system  used  was  a  Sun  SPARC  2  workstation  with 
the  computer  code  written  in  C.  Intermediate  graphics  results  were  created  using  an  in- 
house  computer  graphics  program.  Final  graphics  output  was  performed  by  sending 
output  data  to  the  Deltagraph  graphics  package.  The  sequence  of  steps  of  the  computer 
code  arc  explained  below.  The  computer  code  is  included  as  Appendix  B. 

1.  Initialization 

The  main  computer  program  is  compiled,  along  with  the  header  file  'rkk.h'  and  the 
function  'derivs'  immediately  preceding  it.  The  function  'derivs'  contains  Mingori's 
equations  of  motion  rewritten  into  Equation  (145),  suitable  for  numerical  integration. 
Included  in  the  header  file  'rkk.h'  is  the  numerical  integration  routine  'rk4,'  the  adaptive  step 
size  function  'rkqc,'  and  the  driver  for  the  numerical  integration  routine  'odeint.'  Also 
included  are  functions  for  creating  and  freeing  the  vectors  and  matrices  used  by  the 
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conqjuter  to  store  the  data,  and  a  function  for  error  messages. 

a.  Variables 

The  variables  required  for  both  the  numerical  integration  routine  and 
subsequent  calculations  are  defined  as  global  variables  before  the  main  program.  Variables 
used  only  within  a  specific  function  or  only  in  the  main  program  are  defined  in  their 
respective  functions.  Symbolic  constants  are  also  defined  before  the  main  program  using 
the  #define  statement 

The  time  dependent  variables  of  motion  of  the  Mingori  system, 
CDi,  Q>2,  fiJ3,  2,  i,  z\  z\  are  stored  in  a  seven  column  matrix.  All  other  parameters  to  be 
graphed  are  stored  in  one  dimensional  vectors.  Storage  in  the  vectors  and  the  matrix 
permits  retrieval  by  the  graphics  subroutine  for  plotting  the  variable  over  time. 

b.  Input  File 

The  main  computer  program  scans  the  input  file  for  the  necessary  parameters. 
All  parameters  of  the  Mingori  system,  /i,  /2,  /s,  Af,  m,  fl,  it,  c,  I\\  li,  h\  M\m\  a\  k\ 
c\  and  L,  are  controlled  with  the  input  file.  Also,  the  initial  conditions  for  the  time 
dependent  variables  of  the  system,  coi,  fi>2,  fi>3,  z,  z,  z',  z',  are  specified.  The  length  of 
time  of  the  simulation  and  a  variable  determining  the  desired  accuracy  are  specified. 
Finally,  the  exponential  factor  and  the  final  energy  for  the  postulated  core  energy  functions 
are  read. 

2.  Preliminary  Calculations 

After  the  input  values  are  read  in,  initial  calculations  are  performed  prior  to  the 
numerical  integration  routine.  All  of  the  defined  terms  used  by  D.  L.  Mingori  [Ref  3]  in  his 
non-linear  equations  of  motion  are  computed,  as  well  as  definitions  required  for  the  core 
energy  calculations,  U,  U\  U  totals  U>  h\  U  total,  and  hrigid- 

An  option  is  available  to  specify  critical  damping  in  either  the  platform  or  the 
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rotor.  By  using  the  number  1000.0  in  the  input  file  for  the  damping  coefficient,  the  main 
program  will  automatically  compute  the  coefficient  required  for  critical  damping  of  the 
mass-spring-dashpot  system  and  then  use  this  value  in  all  subsequent  calculations. 

Also  computed  is  the  factor  'dxsaV,  used  in  determining  when  to  save  data.  The 
steps  between  evaluating  the  equations  of  motion  can  become  small,  particularly  when  high 
accuracy  is  desired.  The  interval  required  for  graphics  resolution  is  not  as  restrictive. 
Accordingly,  the  variables  are  saved  only  if  the  step  is  greater  than  the  previously  saved 
step  by  the  factor  'dxsav.' 

3.  Numerical  Integration 

Mingori's  nonlinear  equations  of  motion  arc  numerically  integrated  by  the  fourth 
order  Runge-Kutta  method,  with  adaptive  step  size  control.  The  computer  code  used  is 
based  on  the  Runge-Kutta  method  listed  in  Numerical  Recipes  in  C,  [Ref  8].  The  adaptive 
step  size  control  permits  larger  integration  steps  during  smooth,  well  behaved  portions  of 
the  functions,  and  smaller  steps  during  the  more  irregular  sections  of  the  functions.  The 
integration  routine  accuracy  can  be  controlled  by  a  variable  in  the  input  file. 

Mingori's  non-linear  equations  of  motion  are  contained  in  the  function  'derivs.' 
As  described  previously,  the  equations  have  been  rewritten  as  a  scries  of  first  order  coupled 
differential  equations  suitable  for  numerical  integration. 

The  output  of  the  time  dependent  variables  of  the  system,  fl>i,  0)2,  0)3,  z,  i,  z',  z', 
is  stored  in  an  array  of  seven  columns,  with  the  number  of  rows  required  a  function  of  the 
specified  time  interval  and  accuracy.  A  one  dimensional  vector  is  also  created,  with  the 
time  stored  for  each  step  saved.  This  permits  plotting  the  time  dependent  variables  and 
other  quantities  versus  the  time. 

4.  Calculation  of  System  Parameters 

With  the  array  of  the  time  dependent  variables  of  motion  as  a  function  of  time,  the 
other  system  quantities  listed  in  Section  A  may  now  be  calculated,  with  the  values  stored  in 
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vectors  suitable  for  graphing. 

5.  Graphics  Output 

The  remainder  of  the  program  is  the  necessary  code  for  the  in-house  graphics 
program.  The  graphics  program  is  initialized,  and  the  graphics  window  is  opened.  The 
following  parameters  are  then  plotted  over  time:  Q)i,  0)2.  0)3-0)3  inUiai,  z.  z,  2'.  z', 

»  ~  f^initiaU  ^  ~  ^inUialt  ^  total  ~  '^  total  initials  ^c  ~  ^c  initiatt  ^c  postulated  ~  ^c  postulated  initial 
Ec  —  Ec  initial .  ^c  postulated'  —  Ec  postulated  initial  >  ^»  "Hy  ^^  if.    The  windOW  is  then  cloSCd 

and  the  initial  conditions  and  pertinent  time  dependent  variables  and  other  quantities  of  the 
simulation  are  then  printed.  This  displays  the  graphical  representation  of  the  behavior  of 
the  dual-spin  system,  and  provides  the  actual  initial  and  final  values  of  the  variables  and 
associated  quantities. 

For  final  graphics  output,  the  Deltagraph  graphics  program  is  utilized.  Computer 
simulation  data  is  sent  as  an  output  file.  The  data  is  then  manipulated  into  a  suitable  graphic 
with  proper  scaling  and  axes  limits  to  best  represent  the  dynamics  of  the  computer 
simulation.  The  graphs  are  contained  in  Appendix  A. 

6.  Computer  Program  Validation 

The  two  aspects  of  the  computer  code  requiring  validation  were  the  numerical 
integration  routine  and  Mingori's  equations  of  motion  with  its  associated  dynamical 
quantities  for  the  dual-spin  system. 

The  validation  of  the  numerical  integration  routine  was  performed  by  using  the 
equations  of  motion  for  a  simple  torque-fi^ee  axisymmetric  body.  Various  initial  conditions 
were  read  in,  and  the  dynamical  quantities  were  then  plotted  versus  time.  The  plotted 
behavior  was  then  compared  to  the  actual  response  of  the  system. 

The  validation  of  Mingori's  equations  of  motion  and  the  associated  dynamical 
quantities  required  a  more  comprehensive  approach.  The  equations  of  motion  required 
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validaticMi  for  proper  behavior  of  each  time-dependent  variable  of  motion.  Then  the  kinetic 
energy,  total  energy,  and  angular  momentum  must  be  verified.  Several  cases  were  run 
where  the  platform  mass  and  inertia  would  become  infinitesimally  small.  Then  cases  were 
run  where  the  rotor  mass  and  inertia  would  become  infinitesimally  small.  In  each  of  these 
cases,  subcases  were  run  where  the  mass-spring-dashpot  system  would  be  large, 
dominating  the  dynamics,  to  the  subcase  where  it  became  very  small,  so  the  system 
approximates  a  rigid  body.  These  different  scenarios  would  uncouple  and  isolate  the 
various  parts  of  the  equations  of  motion  to  verify  proper  derivation  of  the  equations.  In 
general,  the  platform  and  the  rotor  were  tested  under  conditions  varying  from  those  in  a 
rigid  body  scenario  to  those  in  a  lightly  damped  body.  The  system  was  then  tested  as  a 
rigid  dual-spin  system,  as  a  dual-spin  system  with  rotor  damping  only,  and  as  a  dual  spin 
system  with  platform  damping  only.  The  dual-spin  system  was  then  tested  with  both  the 
rotor  and  the  platform  containing  dampers.  For  each  case,  the  time  dependent  variables  of 
motion  were  plotted  and  analyzed.  In  all  cases,  the  angular  momentum  was  compared  with 
the  initial  angular  momentum.  Since  all  cases  were  torque  free,  the  angular  momentum 
must  remain  constant.  In  all  cases  explored,  the  angular  momentum  verified  the 
correctness  of  the  equations  of  motion.  The  code  validation  cases  were  not  included  in  this 
thesis  due  to  the  large  number  of  graphs  and  data  that  were  required  to  establish  validation. 
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V.  ANALYSIS 

A.    INTRODUCTION 

1.  Objective 

The  development  of  the  revised  energy-sink  stability  theory  has  been  presented  in 
Chapter  n,  and  the  equations  of  motion  for  a  dual-spin  system  are  contained  in  Qiapter  m. 
With  the  numerical  integration  code  of  Chapter  IV,  the  revised  energy  sink  stability  theory 
can  now  be  verified.  The  core  energy  of  the  system  is  plotted  as  a  function  of  time  and 
compared  to  the  total  energy  for  agreement  with  theory.  The  stability  criterion  computed 
from  the  numerical  simulation  must  then  agree  with  Equation  (97).  Conservation  of 
angular  momentum  is  used  in  each  case  to  verify  the  correcmess  of  the  equations,  and  to 
ensure  the  accuracy  of  the  numerical  integration  routine.  To  approximate  the  core  energies 
over  time,  an  exponential  model  and  a  model  based  on  logistic  growth  are  then  explored. 
The  logistic  growth  model  uses  an  equation  presented  by  Verhulst  [Ref  7],  and  is  referred 
to  as  the  Verhulst  model.  Postulated  models  of  core  energy  and  the  associated  nutation 
angles  are  then  compared  with  the  actual  core  energies  and  nutation  angle  for  agreement 

2.  Numerical  Simulation  Cases 

Four  distinct  cases  needed  to  be  addressed  With  the  inertia  ratio  of  the  dual-spin 
system  greater  than  one,  a  stable  case  and  an  unstable  case  are  analyzed  (Cases  (1)  and 
(4)).  With  the  inertia  ratio  less  than  one,  a  stable  case  and  an  unstable  case  are  also 
analyzed  (Cases  (2)  and  (3)).  For  these  four  cases,  an  exponential  model  is  used  to 
postulate  the  core  energy  and  the  corresponding  nutation  angle.  For  the  case  of  the  inertia 
ratio  greater  than  one  and  unstable,  and  the  two  cases  of  the  inertia  ratio  less  than  one, 
stable  and  unstable,  the  Verhulst  model  is  postulated  (Cases  (5),  (6),  and  (7)).  The  case  of 
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the  inertia  ratio  greater  than  one  and  stable  was  not  included;  from  the  excellent  agreement 
of  Cases  (5),  (6),  and  (7),  one  can  see  that  excellent  agreement  would  also  occur  for  this 
trivial  case,  making  it  unnecessary  to  include.  All  seven  cases  are  summarized  in  Table  (1). 
For  each  case,  the  dynamical  quantities  are  plotted  versus  time  until  the  system  reaches  a 
stable  state.  Each  quantity  is  then  plotted  for  the  first  one  hundred  seconds  to  show  detail. 


h^tal 

^C 

Ed 

CASE 

h  total 

STABILITY 

X 
negative 

X' 

negative 

MODEL 

COMMENTS 

1 

1.039 

stable 

exponential 

large  damping 
in  platform 

2 

0.907 

stable 

negative 

negative 

exponential 

large  damping 
in  platform 

3 

0.897 

unstable 

positive 

positive 

exponential 

insufficient  damping 
in  platform 

4 

1.025 

unstable 

positive 

positive 

exponential 

insufficient  damping 
in  platform 

5 

0.907 

stable 

negative 

negative 

Verhulst 

same  conditions 
asCase2 

6 

0.897 

unstable 

positive 

positive 

Verhulst 

same  conditions 
asCase3 

7 

1.025 

unstable 

positive 

positive 

Verhulst 

same  conditions 
as  Case  4 

Table  (1)  Summary  of  Analyzed  Cases 


In  Appe^  "\  A,  tables  (2)  through  (8)  lis     le  system  parameters,  core  energy 
parameters,  and  conditions.  Immediately  folk    ing  these  tables  are  the  graphs  of  the 

important  dynar.  j.  quantities.  A  typical  geosynchronous  dual-spin  satellite  was  selected 
for  the  numerical  simulation.  The  platform  and  rotor  masses  and  inertias  are  listed  as  part 
of  the  initial  conditions,  and  are  the  same  for  all  cases.  The  inertia  ratio  is  made  greater 
than  or  less  than  one  by  selecting  L,  the  distance  between  the  rotor  and  platform  centers  of 
mass,  to  be  0.3  or  1.0  meter  respectively.   As  a  default  set  of  values,  all  of  the  mass- 

56 


spring-dashpot  systems  have  m  and  m  equal  to  1.0  kg,  k  and  W  equal  to  1.0  N/m,  and  c 
and  c'  equal  to  1.0  kg/sec.  To  establish  the  stable  cases,  additional  energy  dissipation  is 
required  in  the  platform.  To  achieve  this,  cases  (1),  (2),  and  (5)  have  m  increased  to  20.0 
kg  and  c'  increased  to  10.0  kg/sec.  In  all  seven  cases,  the  initial  conditions  are  the  same. 
The  platform  rotates  at  a  geosynchronous  angular  velocity,  the  rotor  spins  at  the  higher  rate 
of  1.5  rad/sec,  and  a  perturbation  is  introduced  by  an  initial  transverse  angular  velocity  of 
0.1  rad/sec.  The  mass-spring-dashpot  systems  have  no  initial  displacement  or  initial 
velocity.  The  core  energy  parameters  are  also  listed  in  the  tables.  The  initial  core  energies 
are  determined  by  the  system's  initial  conditions.  The  final  core  energies  were  taken  from 
the  numerical  simulation  data.  The  exponential  factors,  r  and  /,  were  then  determined 
through  an  iterative  process  to  best  fit  the  modeled  core  energy  to  the  actual  core  energy. 

B.    DISCUSSION 

The  individual  cases  can  now  be  analyzed.  The  graphs  of  the  dynamical  quantities 
are  explored,  and  the  data  will  affirm  the  revised  stability  theory. 

1.  Angular  Momentum 

For  each  case,  the  angular  momentum  is  plotted  versus  time.  Figures  (5),  (6), 
(16),  (17),  (27),  (28),  (38),  (39),  (49),  (50),  (60),  (61),  (71),  and  (72).  Because  the  dual- 
spin  system  has  no  external  forces,  angular  momentum  must  be  conserved.  For  the  stable 
cases,  there  is  excellent  agreement,  with  angular  momentum  varying  by  less  than  one  one- 
hundredth  of  a  percent  over  the  length  of  the  data  run.  This  confirms  the  equations  of 
motion  and  validates  the  accuracy  of  the  numerical  integration  routine.  For  the  unstable 
cases,  the  angular  momentum  percent  difference  increases  to  approximately  four  one- 
hundredths  of  a  percent.  This  is  attributed  to  the  increased  dynamics  of  the  system  as  it 
establishes  the  spin  about  the  transverse  axis,  introducing  very  small  errors  in  the  numerical 
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integration  routine.  The  errors  remain  very  small,  and  the  equations  of  motion  and  the 
numerical  integration  routine  remain  accurate. 

2.  Total  Energy,  Platform  and  Rotor  Core  Energy 

The  graphs  of  energy  are  in  Figures  (7)  through  (10),  (18)  through  (21),  (29) 
through  (32),  (40)  through  (43),  (51)  through  (54),  (62)  tiirough  (65),  and  (73)  through 
(76).  Several  observations  can  be  made. 

The  total  energy  curve  reaches  equilibrium  before  both  the  platform  and  the  rotor 
core  energies  do.  For  Case  2  and  5,  the  total  energy  appears  to  reach  a  maximum  and  then 
drop  down  to  a  final  equilibrium  value.  Careful  comparison  of  Figures  (18)  and  (49), 
shows  the  same  system  with  the  same  initial  conditions,  but  with  a  slightiy  different  curve. 
It  can  be  deduced  that  the  energy  is  not  steady,  but  is  still  oscillating.  The  sampling 
frequency  coincidentally  saved  values  near  the  same  magnitude  of  energy  in  the  region 
from  about  3000  seconds  to  10000  seconds. 

The  total  energy  for  Cases  (1),  (3),  (5),  and  (6)  decreased.  For  Cases  (2),  (4), 
and  (7),  representing  both  stable  and  unstable  systems,  it  increased.  This  can  be  explained 
easiest  with  the  use  of  equations  (26)  and  (27).  For  Cases  (2),  (4),  and  (7),  the  system  is 
settiing  out  about  the  axis  with  the  minimum  moment  of  inertia.  Because  there  are  no 
external  forces,  the  angular  momentum  is  constant  In  the  equation  for  angular  momentum. 
Equation  (26),  the  inertia  terms  are  squared.  But  in  the  equation  for  energy.  Equation  (27), 
the  inertia  terms  are  not  squared.  Therefore,  as  the  angular  velocity  transfers  to  the  axis 
with  the  minimum  moment  of  inertia.  Equations  (26)  and  (27)  show  that  the  energy  will 
increase.  These  equations  apply  to  a  simple  dual-spin  system.  Equations  (123)  and  (136) 
for  the  Mingori  dual-spin  system  would  show  the  same  result,  provided  the  energy 
absorbed  by  the  mass-spring-dashpot  system  and  the  energy  associated  with  the  motor  is 
less  than  the  energy  increase  associated  with  transferring  the  spin  to  the  axis  of  minimum 
moment  of  inertia.  This  is  the  situation  for  Cases  (2),  (4),  and  (7). 
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The  total  energy  value  will  always  be  between  that  of  the  platfonn  core  energy 
and  the  rotor  core  energy.  The  platform  core  energy  will  have  smaller  values  as  its 
equation  assumes  that  the  rotor  is  rotating  at  the  same  rate  as  the  platform,  thereby  not 
accounting  for  the  large  amount  of  kinetic  energy  associated  with  the  rotor.  The  rotor  core 
energy  will  be  higher  than  the  total  energy  as  it  assumes  that  the  slowly  spinning  platform 
is  rotating  at  the  same  rate  as  the  rotor,  providing  the  system  with  additional  kinetic  energy. 

The  hundred-seconds  graphs  of  total  energy  versus  time  shows  curves  with 
several  different  behaviors.  The  time  rate  of  change  of  total  energy  includes  the  energy 
dissipation  rate  of  the  rotor  and  platform  mass-spring-dashpot  systems,  and  the  rate  of 
work  due  to  the  motor  torque  maintaining  the  relative  rotation  rate.  In  the  hundred-seconds 
graphs,  aside  from  the  general  slope  of  the  curve,  there  is  no  apparent  correlation  between 
the  specific  behavior  of  the  total  energy  to  that  of  the  core  energies.  Any  relationship  that 
exists  is  masked  by  the  motor  and  mass-spring-dashpot  system's  influences  on  total 
energy. 

3.  Stability  Criterion 

The  revised  stability  criterion  of  equation  (97)  states  that  the  time  rate  of  change  of 
the  core  energy  over  the  respective  nutation  frequency  must  be  less  than  or  equal  to  zero. 
The  sign  of  the  stability  criterion  for  each  case  must  be  determined.  The  time  rate  of  change 
of  the  core  energy  was  determined  by  dividing  the  final  less  the  initial  value  of  the  core 
energy  by  the  length  of  time  of  the  case.  The  nutation  frequencies  were  computed  using 
initial  conditions.  The  sign  of  the  stability  criterion  is  listed  in  Tables  (1)  through  (8).  The 
numerical  value  was  not  listed  because  the  above  assumptions  used  during  the  computation 
give  it  no  merit.  Cases  with  inertia  ratios  very  close  to  one  were  intentionally  selected  to 
test  the  inertia  ratio  in  this  transition  region.  In  all  cases  analyzed,  the  sign  of  the  revised 
stability  criterion  was  consistent  with  the  stability  of  the  system. 
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4.  Postulated  Core  Energy  and  Nutation  Angle 

An  explicit  relationship  for  core  energy  as  a  function  of  time  does  not  exist  If  an 
equation  describing  core  energy  did  exist,  then  by  using  Equations  (101)  and  (106)  the 
nutation  angle  as  a  function  of  core  energy  and  time  could  be  predicted  for  a  dual- spin 
system.  This  could  then  be  extended  to  an  actual  dual-spin  satellite.  Given  a  sufficient 
model  of  the  satellite,  the  stability  and  the  nutation  angle  as  a  function  of  time  could  be 
predicted.  An  objective  of  this  thesis  is  to  see  if  an  equation  for  the  core  energy  can  be 
developed  that  would  adequately  describe  the  nutation  angle  as  a  function  of  time  for  both 
the  stable  and  unstable  conditions.  The  exponential  model  and  the  Verhulst  logistic  model 
were  explored. 

a.  Exponential  Core  Energy  Model 

Observation  of  the  core  energy  as  a  function  of  time  for  a  stable  system  would 
lead  one  to  conclude  that  it  behaves  in  an  exponential  manner.  Equations  (140)  and  (141) 
are  exponential  representations  of  rotor  and  platform  energy  as  a  function  of  time.  The 
initial  core  energy  was  determined  by  the  initial  conditions  of  the  dual-spin  system.  The 
final  core  energy  was  determined  by  running  the  numerical  simulation  and  calculating  it  at 
the  end  of  the  simulation.  If  this  was  not  available,  a  value  could  be  estimated.  Applying 
the  principle  of  conservation  of  angular  momentum  and  noting  that  the  system  will 
eventually  spin  about  one  of  the  primary  axes,  then  the  equations  for  angular  momentum 
and  total  energy  can  be  used  to  solve  for  the  angular  velocity  about  that  axis.  Substituting 
this  into  Equation  (47)  or  (48),  one  would  arrive  at  an  estimated  fmal  core  energy. 
Although  it  would  not  be  the  actual  fmal  core  energy,  as  motor  torque  contribution  and  the 
mass-spring-dashpot  system's  energy  dissipation  was  not  accounted  for,  it  would  be 
sufficiently  close  to  satisfy  the  computational  requirements.  The  exponential  factor  is 
dependent  upon  the  system  parameters  and  initial  conditions.  For  the  cases  presented, 
different  values  were  tried  until  good  agreement  was  established  with  the  core  energy 
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curve.  Additional  research  would  be  required  to  detennine  a  suitable  exponential  factor  for 
an  actual  satellite,  taking  into  account  the  parameters  and  the  initial  conditions.  Cases  1 
through  4  contain  the  numerical  simulation  data  for  the  exponential  model.  The  actual  and 
postulated  curves  for  core  energy  are  contained  in  Figures  (11)  through  (15),  (22)  through 
(24),  (33)  through  (37),  and  (44)  through  (48).  For  cases  1  and  2,  both  stable,  there  is 
excellent  agreement  between  the  core  energy  and  the  postulated  core  energy.  This  in  turn, 
results  in  excellent  agreement  between  the  actual  nutation  angle  and  the  modeled  nutation 
angle.  The  actual  nutation  angle  is  determined  using  the  angular  momentum  quantities,  as 
shown  iquation  (34).  The  time  dependent  variable  required  to  compute  the  nutation 
angle  is  the  angular  velocity  about  the  spin  axis.  The  modeled  nutation  angle  is  determined 
in  Equations  (101)  and  (106),  and  is  a  function  of  only  one  time-dependent  variable,  core 
energy.  Herein  lies  the  potential  of  the  core  energy  theory.  A  sufficient  model  of  core 
energy  over  time,  as  in  Cases  1  and  2,  will  provide  an  excellent  prediction  of  nutation 
angle,  without  requiring  any  knowledge  of  the  specific  angular  velocities  of  the  system  as  a 
function  of  time. 

Cases  3  and  4  illustrate  the  exponential  model  for  an  unstable  dual-spin 
system.  The  exponential  model  for  core  energy,  and  its  associated  nutation  angle,  rapidly 
approach  their  final  values.  The  actual  core  energy  and  nutation  angle,  however,  behave 
quite  differendy.  The  initial  conditions  have  the  system  near  an  unstable  equilibrium.  The 
system  moves  from  the  unstable  to  the  stable  equilibrium  slowly  at  first.  It  then  increases 
the  rate  at  which  it  approaches  stable  equilibrium,  passes  through  an  inflection  point,  and 
then  approaches  equilibrium  asymptotically.  It  is  clear  that  the  exponential  model 
represents  this  behavior  poorly.  Verhulst's  logistic  equation  was  then  addressed  to 
determine  its  adequacy  in  modeling  the  dual  spin  system. 
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b.  Verhulst  Logistic  Core  Energy  Model 

The  Verhulst  logistic  equation  was  introduced  in  Chapter  HI.  Figure  (4)  and 
the  associated  description  explains  the  charaaeristics  of  the  equation.  The  general  shape  of 
the  curves  in  Figure  (4)  is  very  similar  to  the  stable  and  unstable  cases  of  the  dual- spin 
system.  Case  5  is  the  Verhulst  model  of  the  Case  2  stable  system.  Figures  (55)  through 
(59)  show  that  the  Verhulst  logistic  equation  can  achieve  excellent  agreement  with  core 
energy  and  nutation  angle,  just  as  the  exponential  model  did. 

C^ses  6  and  7  are  the  same  as  (Dases  3  and  4,  except  the  Verhulst  logistic 
equation  is  used  to  model  the  core  energies.  Figures  (66)  through  (70)  and  (77)  through 
(81)  illustrate  the  modeled  and  actual  core  energies.  Although  the  Verhulst  model  for  rotor 
core  energy  had  excellent  agreement,  it  performs  poorly  when  modeling  the  unstable 
system.  The  rotor  core  energy  begins  at  an  initial  value  and  then  decreases  to  an 
equilibrium  value.  Referencing  Figure  (4),  the  Verhulst  logistic  equation  will  model  the 
rotor  core  energy  exponentially. 

To  take  advantage  of  the  Veiiiulst  logistic  equation's  curve  beginning  near  the 
the  unstable  equilibrium  and  its  progression  to  the  asymptotically  stable  equilibrium,  the 
core  energy  must  increase  over  time.  The  platform  core  energy  behaves  in  this  manner  as 
the  nutation  angle  goes  from  a  small  angle  to  ninety  degrees.  Figures  (68)  through  (70) 
and  (79)  through  (81)  show  the  Verhulst  modeled  and  the  actual  platform  core  energies, 
and  the  modeled  and  the  actual  nutation  angles.  The  agreement  between  the  actual  and 
modeled  energies  and  nutation  angles  was  very  good.  The  general  shape  of  curve  was 
consistent,  with  only  a  slight  deviation  in  the  center  of  the  curve,  and  then  a  small  deviation 
as  the  nutation  angle  approaches  the  equilibrium  value  of  ninety  degrees.  With  this 
agreement,  it  is  established  that  a  core  energy  model  exists  that  can  represent  both  stable 
and  unstable  cases.  For  each  of  the  cases,  core  energy  is  plotted  versus  nutation  the  angle. 
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Figures  (15),  (26),  (37),  (48),  (59),  (70),  and  (81).  The  graphics  package  permitted  only  a 
semi-log  plot  of  rotor  core  energy.  The  plot  did  not  provide  any  insight  on  relationships 
for  the  dynamics  of  the  dual-spin  system.  For  both  the  stable  and  unstable  cases,  the  log- 
log  plot  between  platform  core  energy  and  the  nutation  angle  is  linear  from  ninety  degrees 
until  about  two  degrees.  For  nutation  angles  less  than  two  degrees,  the  platform  core 
energy  asymptotically  approaches  the  equilibrium  value.  Case  5  and  6  is  the  same  system 
with  the  same  initial  conditions,  with  the  exception  of  the  increased  damping  system  in  the 
platform  for  C!ase  5.  The  initial  platform  core  energies  and  nutation  angles  are  very  nearly 
the  same  for  each  case,  as  shown  in  Figures  (57),  (58),  (59),  (68),  (69),  and  (70).  One 
can  conclude  by  this  observation  and  the  defmition  of  nutation  angle  as  a  function  of  core 
energy.  Equation  (101),  that  there  exists  a  continuous  platform  core  energy  versus  nutation 
angle  curve.  The  appearance  would  look  similar  to  the  curve  created  by  splicing  Figures 
(59)  and  (70)  together.  By  varying  one  parameter,  for  example  platform  energy  dissipation 
rate,  the  system  would  progress  along  this  curve  and  achieve  equilibrium  with  a  zero 
degree  nutation  angle  or  achieve  equilibrium  with  a  ninety  degree  nutation  angle.  This  is 
what  was  done  with  cases  5  and  6.  Adding  as  the  third  dimension  to  the  curve  the  time  rate 
of  change  of  core  energy  would  then  reveal  the  stability,  the  initial  direction,  and  the  rate  at 
which  the  system  will  arrive  at  the  equilibrium  condition. 

C.   FURTHER  RESEARCH 

The  revised  stability  criterion  for  a  dual-spin,  quasi-rigid,  axisymmetric  system  was 
established.  Numerical  simulation  was  then  used  to  verify  the  revised  stability  theory. 
Further  research  could  be  conducted  in  several  areas.  The  specific  contributions  to  the  total 
energy  could  provide  some  insight  How  much  energy  and  in  what  manner  does  the  motor 
"■  torque  contribute  to  the  total  energy  for  both  the  stable  and  the  unstable  cases?  Also,  by 
plotting  the  energy  dissipation  system's  contributions  over  time,  one  could  determine  its 
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effect  CMi  total  energy,  as  well  as  comparing  it  to  the  core  energy  theory  in  Equation  (97). 

Relationships  for  core  energy  as  a  function  of  time  were  explored.  It  was  established 
that  the  Verhulst  logistic  equation  could  be  applied  to  a  stable  dual-spin  system.  By 
changing  its  parameters,  this  equation  also  applied  to  the  unstable  dual-spin  system  with 
very  good  agreement.  Further  research  could  be  directed  to  find  one  equation  that  would 
address  both  the  stable  and  unstable  cases  of  the  dual-spin  system  without  changing  its 
parameters.  The  equation  for  logistic  growth  with  a  threshold  [Ref.  7]  as  applied  to  the 
platform  core  energy  shows  promise  in  this  regard. 

Ec'  =  V  1 ^^1  1 ^^—\^c  (146) 

V        Ec  unstable  l\        Ec final  I 

Given  the  initial  conditions,  the  platform  core  energy  will  typically  start  at  some 
intermediate  value  and  will  either  increase  or  decrease  as  the  dual-spin  system  reaches 
equilibrium  with  a  nutation  angle  of  either  zero  degrees  or  ninety  degrees.  Equation  (146) 
would  be  well  suited  to  model  the  platform  core  energy.  There  exists  a  platform  core 
energy  value,  Ec  unstable  r  where  the  system  is  at  an  unstable  equilibrium.  Equation  (146) 
shows  that  at  exactly  this  value,  the  rate  of  change  of  the  platform  core  energy  will  be  zero. 
For  any  value  below  the  unstable  equilibrium  value,  the  core  energy  will  jq)proach  the  value 
of  zero.  For  Case  5,  the  final  core  energy  for  the  stable  condition  was  0.122  J.  Although 
this  final  condition  cannot  be  represented  in  Equation  (146),  it  is  sufficiently  close  that  the 
equation  may  still  be  used  to  represent  the  platform  core  energy.  For  any  value  above  the 
unstable  equilibrium  value,  the  core  energy  will  approach  the  equilibrium  value  of  £c/ma/. 
associated  with  the  nutation  angle  at  ninety  degrees.  Initial  conditions  would  determine  the 
initial  core  energy.  The  final  core  energy  can  be  estimated  as  described  earlier  in  this 
chapter.  Finally,  the  exponential  factor  /  may  then  be  determined  based  on  the  system 
parameters  and  initial  conditions. 

Additional  analysis  can  be  performed  on  the  revised  stability  criterion.  The  time  rate 
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of  change  of  core  energy  over  nutadonal  frequency  may  be  calculated  and  plotted  versus 
time  or  versus  other  parameters.  The  behavior  of  this  stability  criterion  and  the  magnitude 
of  it  for  various  conditions  could  provide  some  insight  on  postulating  a  relationship  for 
core  energy  as  a  function  of  time. 


65 


VI.  CONCLUSION 


The  existing  energy-sink  stability  criterion  was  introduced  as 


(41) 


It  was  shown  that  an  inconsistency  in  the  development  disproves  the  existing 
assumption  that  the  motor  energy  input  exactly  balances  shaft  firictional  losses.  A  revised 
energy-sink  stability  criterion  was  then  developed  based  on  Hubert's  definition  of  core 
energy  and  was  presented  as 


(97) 


This  criterion  compliments  the  existing  theory.  Numerical  simulation  was  required  to 
validate  the  theory.  The  Mingori  dual-spin,  quasi-rigid,  axisymmetric  system  was  selected 
for  the  numerical  simulation.  Several  cases  were  analyzed  to  verify  the  revised  energy-sink 
stability  criterion.  By  correctly  postulating  the  platform  or  the  rotor  core  energy,  the 
stability  of  the  system  could  be  determined.  Specific  knowledge  of  the  energy  dissipation 
rates  for  both  the  platform  and  the  rotor  are  no  longer  required. 

An  exponential  model  and  the  Verhulst  logistic  model  for  core  energy,  and  their 
relationship  to  the  nutation  angle,  were  explored.  The  exponential  model  had  excellent 
agreement  with  the  stable  cases,  but  was  inadequate  in  representing  the  unstable  cases.  The 
Verhulst  logistic  model  established  that  an  explicit  relationship  for  core  energy  could  be 
developed.  Nutation  angle  as  a  function  of  core  energy  for  the  dual-spin  system  could  then 
be  predicted.  The  excellent  agreement  of  the  postulated  core  energy  and  nutation  angle  with 
the  actual  core  energy  and  nutation  angle  confirms  the  revised  energy-sink  stability 
criterion.  Additional  research  is  required  to  find  an  optimum  equation  for  core  energy  as  a 
function  of  time  to  represent  both  the  stable  and  unstable  cases. 
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APPENDIX  A    NUMERICAL  SIMULATION  DATA 


CASE  1:  ^  >  1,  Stable,  Exponential  Model 


System 

Parameters                                                   Initial  Conditions 

Platform 

Rotor 

Platform                       Rotor 

/3'  = 

1600  kgm2 
1500  kgm2 

1000  kgm2 
1200  kgm2 

z'  =          0.0  m              z    =   0.0  m 

z'  =          0.0  SI              i    =   0.0  SI 
s                              s 

M'  = 

1000.0  kg 

M  = 

700.0  kg 

G>3'=  7.27x10-5^       fl>j=  1.5^ 

m  = 

20.0  kg 
1.0  m 

m  = 
a  = 

1.0  kg 
1.0  m 

coi=       0.10^            ^=00^ 

r  = 

10  S 
m 

k  = 

1-0  S 
m 

Core  Energy  Parameters 

c'  = 

10.0  j'i 

sec 

c  = 

1.0  ^^ 

sec 

Platform                            Rotor 

Ecinitia{=         13.402  J               EcinUial^         3145.4  J 

L  = 

0.3  m 

h  total 
h  total 

=      1.039 

Ecfinal      =          0.072  J                Ec  final     =         3161.7  J 

Ec                                        Ec 

—     =        negative             — ^     =        negaave 

/      =    -.00134  s-i              r      =    -.00134  s- 

Table  (2)   Case  1  Parameters 
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Case  1:  Is/It>l,  Stable,  Exponential  Model 
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Figure  (5)   Total  Energy  and  Percent  Difference  Angular  Momentum  Versus  Time 
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Case  1:  Is/It>l,  Stable,  Exponential  Model 
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Figure  (6)   Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  100  Seconds 
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Case  1:  Is/It>l, Stable,  Exponential  Model 
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Figure  (7)   Total  Energy  and  Rotor  Core  Energy  Versus  Time 
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Case  1:  Is/It>l, Stable,  Exponential  Model 
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Figure  (8)   Total  Energy  and  Rotor  Core  Energy  First  100  Seconds 
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Case  l:Is/It>l,  Stable,  Exponential  Model                             | 
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Figure  ( 1 0)   Total  Energy  and  Platform  Core  Energy  -  First  100  Seconds 


71 


3162 


Case  l:Is/It>l,  Stable,  Exponential  Model 
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Figure  (11)   Modeled  and  Actual  Rotor  Core  Eno-gy  and  Nutation  Angle  Vctsus  Time 
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Case  1:  Is/It>l,  Stable,  Exponential  Model 
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Figure  ( 12)  Modeled  and  Actual  Rotor  Core  Energy  and  Nutation  Angle-  First  100  Seconds 
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Case  1:  Is/It>l,  Stable,  Exponential  Model 
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Figure  ( 13)   Modeled  and  Actual  Platform  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  1:  Is/I t> Instable,  Exponential  Model 
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Figure  ( 14)  Modeled,  Actual  Platform  Core  Energy  and  Nutation  Angle-  First  100  Seconds 
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Case  1:  IsAt>l,  Stable,  Exponential  Model 
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Figure  (15)    Core  Energy  Versus  Nutation  Angle 
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CASE  2:  ^  <  1  ,  Stable,  Exponential  Model 


System 

Parameters                                                 Initial  Condiocos 

Platform 

Rotor 

Platform                      Rotor 

f3  = 

1600  kgm2 
1500  kgm2 

/3  = 

1000  kgm2 
1200  kgm2 

z'  =          0.0  m               z    =    0.0  m 

z'  =         0.0  n^            z  =  0.0  ni 
s                              s 

lf'  = 

1000.0  kg 

A/  = 

700.0  kg 

fi)3'=  7.27  X  10-5  mi      6,3  =  1.5  mi 

sec          ^           sec 

/ 
71   = 

20.0  kg 
1.0  m 

m  = 

1.0  kg 
1.0  m 

o>i  =        0.10  5^             G)2=  0.0^ 

sec               ^           sec 

fl'  = 

fe'  = 

10  S 
m 

ifc  = 

1-0  S 
m 

Core  Energy  Parameters 

c'  = 

10.0^ 

c  = 

10^ 

Platfonii                            Rotor 

sec 

sec 

Ec  initial' =           15.34  J                Ecimtial=         3147.3  J 

L  = 

1.0  m 

'j  ro/a/ 

=     0.907 

Ec  final     =          0.123  J               Ecfiruil    =         3170.9  J 

Ec                                       Ec 

=        negative             — ^     =         negative 

X'                                       X 

/      =     -.00102  s-i              r      =     -.00102  s-i 

Table  (3)   Case  2  Parameters 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  ( 1 6)   Total  Energy  and  Percent  Difference  Angular  Momentum  Vctsus  Time 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (17)   Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  100  Seconds 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (19)   Total  Energy  and  Rotor  Core  Energy  -  First  100  Seconds 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (20)   Total  Energy  and  Platform  Core  Energy  Vctsus  Time 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (21)    Total  Energy  and  Platform  Core  Energy  -  First  100  Seconds 
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3180 


Case  2:Is/It<l,  Stable,  Exponential  Model 
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Figure  (22)   Modeled  and  Actual  Rotor  Core  EnCTgy  and  Nutation  Angle  Versus  Time 
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Case  2;  Is/It<l,  Stable,  Exponential  Model 
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Figure  (23)    Modeled,  Actual  Rotor  Core  Energy  and  Nutation  Angle  -  First  100  Seconds 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (24)    Modeled  and  Actual  Platform  Core  EnCTgy  and  Nutation  Angle  Versus  Time 


16 


\5S 


3   15 

00 

w 

c 

^  ]45 


14 


13.5 


Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (25)  Modeled,  Actual  Platform  Core  Energy  and  Nutation  Angle-  First  100  Seconds 
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Case  2:  Is/It<l,  Stable,  Exponential  Model 
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Figure  (26)   Core  Energy  Versus  Nutation  Angle 
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CASE  3:  ^  <  1  ,  Unstable,  Exponential  Model 


System  Parameters 
Platform  Rotor 


Ii=  1600  kgm2 
l3  =  1500  kgm2 
A/'=  1000.0  kg 


m  = 
fl'  = 
k'  = 

c  = 


1.0  kg 
1.0  m 


10  S- 
m 


1.0 


N 
m 

sec 


h  = 

M  = 
m  = 
a- 
k  = 

c  = 


1000  kgm2 

1200  kgm2 

700.0  kg 

1.0  kg 

1.0  m 


10  S 
m 


1.0 


N 
m 

kg 
sec 


L=       1.0  m         lnmL=     0.897 

h  total 


Initial  Conditions 
Platform  Rotor 


/  = 

0.0  m 

z   = 

0.0  m 

z'  = 

0.0  m 

z    = 

0.0  ni 

S 

6)3- = 

7.27  X  10-5  rad 
sec 

0)3  = 

1.5^^ 
sec 

0)1    = 

0.10^ 
sec 

0)2  = 

0.0^^ 
sec 

Core  Energy  Parameters 
Platform  Rotor 


^c  initial  ~ 

r  '     _ 

^cfinal      ~ 

Ec_ 
X 


15.1  J 
1160.9  J 

positive 
-.0005  s-i 


^c  initial  ~ 
^c  final    ~ 

Ec_ 

A 

r      = 


3061.6  J 
1489.9  J 

positive 
-.0005  s-i 


Table  (4)   Case  3  Parameters 
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Case  3:  Is/lKl,  Unstable,     Exponential  Model 
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Figure  (27)   Total  Energy  and  Percent  Difference  Angular  Momentum  Versus  Time 
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Case  3:Is/It<l,  Unstable,  Exponential  Model 
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Figure  (28)    Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  1(X)  Seconds 


83 


1370 

Case  3: Is/It<l,  Unstable,  Exponential  Model 
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Figure  (29)   Total  Energy  and  Rotor  Core  Energy  Versus  Time 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (30)   Total  Energy  and  Rotor  Core  Energy  -  First  ICX)  Seconds 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (3 1 )   Total  Energy  and  Platform  Core  Energy  Vctsus  Time 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (32)    Total  Energy  and  Platform  Core  Energy  -  First  100  Seconds 
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Case  3:  Is/It<l,  Unstable,  Exponential  Mcxlel 
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Figure  (33)   Modeled  and  Actual  Rotor  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (34)    Modeled,  Actual  Rotor  Core  Biergy  and  Nutation  Angle  -  First  100  Seconds 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (35)   Modeled  and  Actual  Platform  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (36)  Modeled,  Actual  Platform  Core  Energy  and  Nutation  Angle-  First  100  Seconds 
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Case  3:  Is/It<l,  Unstable,  Exponential  Model 
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Figure  (37)    Core  Energy  Versus  Nutation  Angle 
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CASE  4:  ^  >  1   ,  Unstable,  Exponential  Model 


System  Parameters 
Platform  Rotor 


Initial  Conditions 
Platform  Rotor 


/l'= 

1600  kgm2 

/i  = 

1000  kgm2 

/3'  = 

1500  kgm2 

/3  = 

1200  kgm2 

A/'  = 

1000.0  kg 

M  = 

700.0  kg 

m  = 

1.0  kg 

m  = 

1.0  kg 

a'  = 

1.0  m 

a  = 

1.0  m 

k'  = 

10  S 
m 

k  = 

1-0  S 
m 

c'  = 

1.0  ^^ 

sec 

c- 

1.0  ^^ 

sec 

L  = 

0.3  m 

h  total  , 
It  total 

=      1.025 

z'  = 

0.0  m 

z   = 

0.0  m 

z'  = 

0.0  ni 

z   = 

0.0  m 

0)3- = 

7.27  X  10-5  ^Sd 
sec 

a>3  = 

^•^^ 

(Ox  = 

0.10^ 
sec 

G)2  = 

0-0  i^ 

Core  Energy  Parameters 
Platform  Rotor 


^c  initial  ~ 

E'  _ 

c  final  ~ 

X 


13.2  J 
1247.1  J 

positive 
-.0005  s-i 


t'c  initial  ~ 

^c  final    ~ 

k.       = 

A 

r      = 


3059.7  J 
1552.5  J 

positive 
-.0005  s-i 


Table  (5)   Case  4  Parameters 
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Case  4:  Is/Iol,  Unstable,  Exponential  Model 
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Figure  (38)   Total  Energy  and  Percent  Difference  Angular  Momentum  Versus  Time 
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Figure  (39)   Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  100  Seconds 
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Case  4:  Is/It>l, Unstable,  Exponential  Model 
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Figure  (40)    Total  Energy  and  Rotor  Core  Energy  Versus  Time 
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Case  4:  Is/It>l,  Unstable,  Exponential  Mode 
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Figure  (41)    Total  Energy  and  Rotor  Core  Energy  -  First  1(X)  Seconds 
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Figure  (42)   Total  Energy  and  Platform  Core  Energy  Versus  Time 
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Case  4:  Is/It>l,  Unstable,  Exponential  Model 


- 

•• 

's. 

Ill      1      1     \ 

1         1 

_/"'                      —       total  energy           1 

-^ 

platform  core  energy 

1              1             1             1              1             1             1 

15.2 

15 

14.8 

14.6 

14.4  S 

14.2  \ 

14   I 

13.8 

13.6 

13.4 

13.2 


10    20    30    40    50    60 

time  (seconds) 


70 


80 


90    100 


Figure  (43)    Total  Energy  and  Platform  Core  Energy  -  First  100  Seconds 
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Case  4:Is/It>l,  Unstable,  Exponential  Model 


Figure  (44)   Modeled  and  Actual  Rotor  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  4:Is/It>l,  Unstable,  Exponential  Model 
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Figure  (45)   Modeled,  Actual  Rotor  Core  Energy  and  Nutation  Angle  -  First  1(X)  Seconds 
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Case  4:  Is/It>l,  Unstable,  Exponential  Model 
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Figure  (46)    Modeled  and  Actual  Platform  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  4:  Is/It>l,  Unstable,  Exponential  Model 
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Figure  (47)  Modeled,  Actual  Platform  Core  Energy  and  Nutation  Angle-  First  100  Seconds 


94 


10000 


c 

e 

o 


1000 


Case  4:  Is/It>l,  Unstable,  Exponential  Model 
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Figure  (48)    Core  Energy  Versus  Nutation  Angle 
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CASE  5:  ^  <  1  ,  Stable,  Verhulst  Model 

It 


System 

Parameters                                              Initial  Conditions 

Platform 

Rotor 

Platform                      Rotor 

1600  kgm2 
1500  kgm2 

h  = 

1000  kgm2 
1200  kgm2 

z'  =          0.0  m               z   =    0.0  m 

z'  =         0.0  m^            z  =  0.0  ni 
s                              s 

M'  = 

1000.0  kg 

M  = 

700.0  kg 

0)3'=  7.27x10-5^       a>3=  1.5^ 

m'  = 
a'  = 

20.0  kg 
1.0  m 

m  = 
a  = 

1.0  kg 
1.0  m 

(Ox=       0.10^            a>2=o.O^ 

k'  = 

10  S 
m 

k  = 

10  S 
m 

Core  Energy  Parameters 

c'  = 

10.0  ^ 

c  = 

••0^ 

Platfonn                            Rotor 

sec 

sec 

Ecinitial' =       15.341  J             Ecinitial=       3147.3  J 

L  = 

1.0  m 

h  total 
h  total 

=     0.907 

Ec  final      =       0.122  J              Ecf^nal    =      3170.9  J 

Ec                                    Ec 

=      negative           —     =      negative 

X'                                    X 
/      =   -.0002  s-i            r      =   -.0002  s-i 

Table  (6)   Case  5  Parameters 
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Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (49)   Total  Energy  and  Percent  Difference  Angular  Momentum  Versus  Time 
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Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (50)   Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  100  Seconds 
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Case  5: Is/It< Instable,  Verhulst  Model 
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Case  5: Is/It<l,  Stable,  Verhulst  Model 
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Figure  (52)    Total  Energy  and  Rotor  Core  Energy  -  First  100  Seconds 
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Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (53)   Total  Energy  and  Platform  Core  Energy  Versus  Time 
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Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (54)    Total  Energy  and  Platform  Core  Energy  -  First  100  Seconds 


99 


3175 


Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (55)   Modeled  and  Actual  Rotor  Core  Energy  and  Nutation  Angle  Versus  Time 


Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (56)   Modeled,  Actual  Rotor  Core  Energy  and  Nutation  Angle  -  First  100  Seconds 
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Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (57)   Modeled  and  Actual  Platform  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  5:  Is/It<l,  Stable,  Verhulst  Model 
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Figure  (58)  Modeled,  Actual  Platfonn  Core  Energy  and  Nutation  Angle-  First  100  Seconds 
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3175 

Case  5:  Is/It<l,  Stable, 

Verhulst  Model 
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Figure  (59)    Core  Energy  Versus  Nutation  Angle 
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CASE  6:  ^  <  1  ,  Unstable,  Verhulst  Model 
I, 


System  Parameters 
Platform  Rotor 


'l'= 

1600  kgm2 

/i  = 

1000  kgm2 

'3'= 

1500  kgm2 

/3  = 

1200  kgm2 

4'  = 

1000.0  kg 

A/  = 

700.0  kg 

n  = 

1.0  kg 

/n  = 

1.0  kg 

3'  = 

1.0  m 

a  = 

1.0  m 

t'  = 

10  S 
m 

k  = 

10  S 
m 

r'  = 

1.0  ''8 

sec 

c  = 

1.0  ''^ 

sec 

L  = 

1.0  m 

=     0.897 

Initial  Conditions 
Platform  Rotor 


z'  =          0.0  m 

Z     = 

0.0  m 

z'  =         0.0  ni 

z    = 

0.0^ 

C03'  =  7.27  X  10-5  gi 
sec 

fi>}  = 

^•^1^ 

0,,  =       0.10^ 

a>2  = 

0.0  i^ 

Core  Energy  Parameters 
Platform  Rotor 


^c  initial  ~ 

E'  _ 

c  final  ~ 

X 


15.1  J 
1161.0J 

positve 
-1 X  10-^  s-i 


^c  initial  ~ 
^c  final    ~ 

Ec_     ^ 

X 

r      = 


3061.6  J 
1489.9  J 

positive 
-IxlO-^s-i 


Table  (7)  Case  6  Parameters 
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Case  6:  Is/It<l, Unstable,  Verhulst  Model 
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Figure  (60)   Total  Energy  and  Percent  Difference  Angular  Momentum  Versus  Time 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (61 )   Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  1(X)  Seconds 
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Case  6;  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (62)   Total  Energy  and  Rotor  Core  Energy  Versus  Time 
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Case  6:  Is/I-      Unstable,  Verhulst  Model 
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Figure  (63)   Total  Energy  and  Rotor  Core  Energy  -  First  100  Seconds 
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Case  6:  Is/It<l, Unstable,  Verhulst  Model 
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Figure  (64)   Total  Energy  and  Platform  Core  Energy  Versus  Time 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (65)    Total  Energy  and  Platform  Core  Energy  -  First  1(X)  Seconds 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (66)   Modeled  and  Actual  Rotor  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (67)   Modeled,  Actual  Rotor  Core  Energy  and  Nutation  Angle  -  First  100  Seconds 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (68)    Modeled  and  Actual  Platfonn  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (69)  Modeled,  Actual  Platform  Core  Energy  and  Nutation  Angle-  First  100  Seconds 
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Case  6:  Is/It<l,  Unstable,  Verhulst  Model 
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Figure  (70)    Core  Energy  Versus  Nutation  Angle 


109 


CASE  7:  ^  >  1   ,  Unstable,  Verhulst  Model 


System 

Parameters                                                      Initial  Conditions 

Platform 

Rotor 

Platform                      Rotor 

/l'= 

1600  kgm2 
1500  kgm2 

/i  = 
h  = 

1000  kgm2 
1200  kgm2 

z'  =          0.0  m              z   =   0.0  m 

z'  =          0.0  m^             z   =   0.0  ^ 
s                               s 

M'  = 

1000.0  kg 

M  = 

700.0  kg 

a>3'=  7.27x10-5^       G>3=  1.5^ 

m  = 

1.0  kg 
1.0  m 

m  = 
a  = 

1.0  kg 
1.0  m 

(Ox  =        0.10^             G)2=  O.OgJ 

^                    sec                ^            sec 

k'  = 

10  S 
m 

k  = 

10  S 
m 

Core  Energy  Parameters 

c'  = 

10  S 

c  = 

10^2 

Platform                           Rotor 

sec 

sec 

Ec  initial   =               13.2  J                      Ecimtial=             3059.7. 

L  = 

0.3  m 

h  total 
h  total 

=      1.025 

Ec  final'     =             1247.1  J                   Ecfinal    =             1552.5. 

Ec                                            Ec 

—     =           positve               — ^     =           positive 

A'                                             X 
r       =    -l.lxlO^s-i              r      =    -1.1x10-^ 

Table  (8)   Case  7  Parameters 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 
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Figure  (71)   Total  Energy  and  Percent  Difference  Angular  Momentum  Versus  Time 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 
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Figure  (72)   Total  Energy  and  Percent  Difference  Angular  Momentum  -  First  100  Seconds 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 
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Figure  (73)   Total  Energy  and  Rotor  Core  Energy  Versus  Time 


}'if\A 
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Figure  (74)    Total  Energy  and  Rotor  Core  Eno-gy  -  First  100  Seconds 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 
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Figure  (75)   Total  Energy  and  Platform  Core  Energy  Versus  Time 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 
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Figure  (76)    Total  Energy  and  Platform  Core  Energy  -  First  100  Seconds 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 
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Figure  (77)   Modeled  and  Actual  Rotor  Core  Energy  and  Nutation  Angle  Versus  Time 
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Figure  (78)    Modeled,  Actual  Rotor  Core  Energy  and  Nutation  Angle  -  First  100  Seconds 
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Case  7:  IsAt>l, Unstable,  Verhulst  Model 
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Figure  (79)   Modeled  and  Actual  Platform  Core  Energy  and  Nutation  Angle  Versus  Time 
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Case  7:  Is/It>l, Unstable,  Verhulst  Model 
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Figure  (80)  Modeled,  Actual  Platform  Core  Energy  and  Nutation  Angle-  First  1(X)  Seconds 
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Case  7:  Is/It>l,  Unstable,  Verhulst  Model 


1400 


3000 

- 

2800 

^ 

modeled  rotor  core  energy 

2600 

^2400 

1 

2  2200 

8 

- 

2000 

- 

1800 

—  platform  core  energy 

1600 

— 

—  modeled  platform  core  energy 

10 
outatioQ  angle  (degrees) 


10000 


1000 


-    100 


10 


100 


Figure  (81)    Core  Energy  Versus  Nutation  Angle 
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APPENDIX  B     COMPUTER  PROGRAM  CODE 


declarations 


•include  <math.h> 
llnclude  <malloc.h> 
linclude  <stdio.h> 

♦define  MAXSTF  1500C0 
♦define  TINY  l.Oe-40 


/  *    orio  i  Ml 
/  *    otlri  ni 


♦define  PGROW  -C.20  /*  krqc  */ 

♦define  PSHRNK  -0.25  /*  krqc  */ 

♦define  FCOR  0.06665666666666666666666666666667  /•  1.0/15.0  krqc*/ 

♦define  SAFETY  0.9  /*  krqc  */ 

♦define  ERRCON  6.0e-4  /*  krqc  */ 


float  **y-0,  'xx-O: 

int  kmax-O,  kount-C; 

float  'xp-O,  ••yp-O,  dxsav-0; 


/*  defining  declaration  rkdumh  •/ 
/*  defining  declaration  odeint  */ 
/*  defining  declaration  odeini  */ 


error  function  *•***♦•*****•****•*•**/ 


void  nrerror  (error_text) 

char  error_text 1 ]  ; 

{ 

void  exit  {) ; 

fprintf  (stderr,  "Nur-erical  Recipes  run-time  error  ...  \n")  ; 

fprintf  (stderr,  "*s\r."',  error_text)  ; 

f  printf  (stderr ,  "  .  .  .  .-.ow  exiting  to  system.  ..  \n")  ; 

exit(l); 

1 


/ ....*»........ 

float  "vectorr (nl, nh) 

int  nl,  nh; 

( 

float  *v; 


vector  function  *****♦*♦••*****♦»*•**' 


V  -  (float  *)malloc ( (unsigned)  (nh-nl+1) *sizeof (float )) ; 
if  (!v)  nrerror  ("allccation  failure  in  vector  ()"); 
return  v-nl; 
) 

/•••••.*••••***..*.*.  matrix  function  ****••**♦•••♦*•*»••  *•♦♦••♦•'• ,' 

float  '"matrix  (r.rl,  r.rh,  ncl,  nch) 

int  nrl,  nrh,  ncl,  nch; 

1 

int  i; 

float  ••.-«; 

n;  -  (float  **)  r.allcc  !  (unsigned)  (nrh-nrl  +  1 )  *sizeof  (float *))  ; 

if  (!m)  nrerror ("aliccation  failure  1  in  matrix ()"); 

m  —  nrl; 

for  (i-nrl; i<-nrh;i+-) 

{ 

m[i]  -  (float  *)  malloc  (  (unsigned)  (nch-ncl  +  1)  *si7.eof  (f  l^il ) 

if  (!m|i;)  nrerror  ("allocation  failure  2  in  matrixO"); 

mii)  —  ncl; 

} 
return  m; 
i 


/' 


free  vector  function  ***************** 


void  free_vector (V, ni, nh) 
float  "v; 
int  nl,  nh; 
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free ((char*)  (v+nl) ) ; 
) 

/•••••**••♦**•*♦•«*•**•*•  free  matrix  function  * 

void  f ree_matrix (m,  nrl,  nrh,  ncl,  nch) 

float  **m; 

int  nrl,  nrh,  ncl,  nch; 

( 

int  i: 

for (i-nrh;i>-nrl;i--)   free  ((char*)  (m(i] +ncl) )  ; 

free ((char*)  (m+nrl) ) ; 

) 


/•••*•*•*•••**••*****••  xkA    function  ***•********•*•****•♦*/ 
void  rk4 ty, dydx, n, x, h, yout, derivs) 

float  y(],  dydx[],  x,  h,  yout [ ) ; 
void  (*derivs)  ()  ; 
int  n; 

{ 

int  i; 

float  xh,  hh,  h6,  "dym,  *dyn,  *dyt,  *yt,  *vectorr(); 

void  f ree_vector ()  ; 

dym  -  vectorr (1, n)  ; 

dyn  -  vectorr (1, n) ; 

dyt  -  vectorr (1, n) ; 

yt  -  vectorr (1, n) ; 

hh  -  h*0.5; 

h6  -  h/6.0; 

xh  -  x+hh; 

for  (i-l;i<-n;i  +  +)   yt[i)  -  y  (i) +hh*dydx  I  i)  ;     /*  first  st.ep  *  / 

(•derivs) (xh, yt, dyt, dydx); 

for  (i-l;i<-n;i++)   yt(i)  -  y [i] +hh*dyt (i) ; 

(•derivs)  (xh,  yt,dyir,,  dyt)  ; 
for  (i-1; i<-n;i++) 
{ 

yt(i)  -  yli]*h*dym(il; 
dyn(i)  -  dyt [ i) +dym[i) ; 
} 

(*derivs) (x+h, yt, dyt, dym) ; 

for  (i-l;i<-n;i++)   yout[i]  -  y (i) +h6* (dydx ( i ) +dyt 1 i] +2 . 0 'dyn 1 i i ) 

f ree_vector (yt, 1, n) ; 

f ree_vector (dyt, 1, n) ; 

f ree_vector (dym, 1 , n)  ; 

f ree_vector (dyn, 1, n) ; 

) 


/•••**•**•***•***<*«..*  r)tdumb  function 


void  r)<dup±)  (vstart,r.var,xl,x2,nstep,  derivs) 

int  nvar, nstep; 

float  vstart [ ) , xl, x2; 

void  (*derivs)  ( ) ; 

{ 

int  i,  k; 

float  X,  h; 

float  "v,  *vout,  *dv,  *vectorr(); 

void  r)(4(),  nrerrorO,  f ree_vector  ( )  ; 

V  -  vectorr (1, nvar ) ; 
vout  -  vectorr  (1, nvar ) ; 
dv  -  vectorr (1, nvar ) ; 
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for  (i-l;i<-nvar; i+*) 

{ 

v(i)  -  vstart (i]; 

yli) [IJ  -  v|i); 

) 
XX [1]  -  xl; 
X  -  xl; 

h  -  (x2-xl) /nstep; 
for  (k-l;k<-nstep;k++) 

'.   { 

(*derivs)  (x, v, dv) ; 
rk4 (V, dv, nvar, x, h, vout,derivs) ; 

if  (x+h  --  x)  nrerror  ("Step  size  too  small  in  routino  RKnt'rin") 
X  +-  h; 
xx(k+l)  -  x; 
for  (i-1; i<-nvar; i++) 
( 

V (i ]  -  vout (i] ; 
yli] [k+1]  -  v(il; 
) 
) 
f ree_vector (dv, 1, nvar) ; 
f ree_vector (vout, 1, nvar)  ; 
f ree_vector (v, 1, nvar)  ; 
) 


/••••••*•*•*••**•**•*•  rkqc  function  ***************** 

void  rkqc(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs) 


/ 


float  y[),  dydx [ ) ,  *x,  htry,  eps,  yscal(),  *hdid,  *hnext; 
void  ('derivs)  {) ; 
int  n; 

{ 

int  i; 

float  xsav,  hh,    h,  temp,  errmax; 

float  'dysav,  "ysav,  *ytemp,  *vectorr{); 

void  rk4(),  nrerror  (),  f ree_vector () ; 

dysav  -  vectorr(l,n); 
ysav  -  vectorr (1, n) ; 
ytemp  -  vectorr  (1, n) ; 
xsav  -  ('x)  ; 

for  (i-1;  i<-n;  i-t-+) 

( 

ysavli]  -  y[i]; 

dysav  [i]  -  dydx (i) ; 

) 
h  -  htry; 
for  (  ;  ;  ) 

( 

hh  -  0.5*h; 

rk4 (ysav, dysav, n,  xsav,  hh,  ytemp,  derivs)  ; 

•x  -  xsav+hh; 

(•derivs)  (*x, ytemp, dydx) ; 

rk4 (ytemp, dydx,  n,  'x,  hh,  y,  derivs) ; 

*x  -  xsav+h; 

if  (*x  --  xsav)  nrerror  ("Step  size  too  small  in  routine  HKO'J") 

rk4 (ysav, dysav, n,  xsav, h, ytemp, derivs)  ; 

errmax  -  0.0; 

for  (i-1; i<-n; i++) 
( 

ytemp li)  -  y (i] -ytemp [i] ; 
temp  -  fabs (ytemp (i ) /yscal (i) ) ; 
if  (errmax  <  temp)   errmax  -  temp; 
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) 

errmax  /-  eps; 

if  (errmax  <-  1.0) 

I 

*hdid  -  h; 

♦hnext  -  (errmax  >  ERRCON  ?  SAFETYMi'oxpCrCPOW  l.-n  (n  i  m  ,- m 

-1  I! 

break ; 

) 
h  -  SAFETyh*exp(PSHRNK*log  (errmax)  )  ; 

) 
for  (i-l;i<-n;i++)   y(l]  +-  ytemp[i ) *FCOR; 
f ree_vector  (ytemp, 1, n)  ; 
f ree_vector (dysav, i, n)  ; 
f ree_vector (ysav, 1, n) ; 
}  _ 

/• ...♦.«.. odeint  function  ..«.♦•♦*••*•»»***♦»••/ 

void  odeint  (y start,  nvar,xl,  x2,eps,  hi,  hmin,  no)c,  nbad,derivs,  rVqi-) 

float  ystarti],  xl,  x2,  eps,  hi,  hmin; 
int  nvar,  *no)c,  *nbad; 
void (*derivs) () ; 
void(*r)^qc)  ()  ; 

I 

int  nstp,  i; 

float  xsav,  x,  hnext,  hdid,  h; 
float  •yscal,  "y,  'dydx,  *vectorr () ; 
void  nrerrorO,  f  ree_vector  ()  ; 

yscal  -  vectorr  (1, nvar) ; 
y  -  vectorr  (1, nvar) ; 
dydx  -  vectorr (1, nvar)  ; 
X  -  xl; 

h  -  (x2  >  xl)  ?  fabs{hl)  :  -fabs(hl); 
•no)c  -  (*nbad)  -  Jco'jnt  -  0; 
for  (i-1;  i<-nvar;  i+-')   y(i]  -  ystart(i); 
if  ()cmax  >  0)   xsav  -  x-dxsav*2.0; 
for  {nstp-l;nstp<-MAXSTP;nstp++) 
( 

(*derivs)  (x, y, dydx) ; 

for  (i-1;  i<-r.var;  i  +  +  )   yscalli]  -  f  abs  (y  I  i  )  )  +  f  abs  (dydx  I  i  ]  '  h  )  •  7  :  :JY; 
if  (kmax  >  C) 
( 

if  ;f  abs  (x-xsav)  >  f  abs  (d:<sav)  ) 
{ 

if  (kount  <  kmax-1) 
{ 

xp(++kount)  -  x; 

for  (i-1; i<-nvar; i  +  +  )   yp  1  i )( kount ] -y i i ] ; 
xsav  -  x; 
) 
) 
) 
if  (  (x+h-x2/ ' (x+h-xl)  >  0.0)   h  -  x2-x; 

( *rkqc)  (y, dycx, nvar, tx, h, eps, yscal, shdid,  shnext, dcrivs) ; 
if  (hdid  ~  h)  ++(*nok);   else  ++ Cnbad) ; 
if  (  (x-x2)*  (y.2-xl)  >-  0.0) 
{ 
for  'i-1; i<-nvar; i++) 

ystart(i)  -  y [i] ; 
if  (kmax) 
( 

xp(++kount )-x; 
for  (i-1;  i<-nvar;  i  +  +  )   yp[i]  (kount)  "■    yMI; 
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) 

f ree_vector (dydx, 1, nvar) ; 
f ree_vector (y, 1, nvar) ; 
f ree_vector (yscal, 1, nvar) ; 
return; 
) 
if  (fabs(hnext)  <-  hmin)  nrerror  ("Step  size  too  small  in  oni:irii") 
h  -  hnext; 
) 
nrerror ("Too  many  steps  in  routine  ODEINT") ; 
) 
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/*•**•••**«•«*•*•••••♦  declarations  ••*•*•**•*****♦  »«.««..♦..  •«  •  •  .  •  ' , 
linclude  "rkk.h" 

♦define  N  7 

♦define  MAXARRAY  100000 

•  define  SQ(x)  ((x)*  (x) ) 

int  kount,  kmax; 

float  **yp,  *xp; 

float  dxsav,  dxmin,  dxllmit; 

float  Il-O.O,  12-0.0,  13-0.0,  M-0.0,  m-0.0,  mb-0.0,  a=0.0,  k-O.O,  r-n.'>; 

float  Ilp-0.0,  I2p-0.0,I3p-0.0,Mp-0.0,mp-O.0,mbp-0.0,ap-0.0,kp-0.o,'-,.  no. 

float  L-0.0,  sigma-O.O,  HT-0.0,  v-0.0,  rho-0.0,  rhop-0.0,  Ll-O.i),  l.?-(i.(): 

float  A-0.0,  C-0.0,  J3p-0.0,  w3p-0.0; 


/••••***•«***•***•******.  EQUATIONS  OF  MOTION  *********** , 

void  derivs (x, y,dydx) 
float  X,  y[],  dydxU; 
( 

float  zeta-G.O,  dydxzeta-0 .0; 

float  Al-O.O,    A2-0.0,   A3-0.0,  A4-0.0,   A5-0.0,   A6=0.0,   A7-0.0,   AP-O.O. 

float  A9-0.0,   AlO-O.O,  All-0.0,  A12-0.0,  A13-0.0,  A11=0.0,  A13-0.0,  Al^-n.o. 

float  A17-0.0,  A18-0.0,  A19-0.0,  A20-0.0,  A21-0.0,  A22-0.0,  A23-n.O,  A^l-n  0. 

float  A25-0.0,  A26-0.0,  A27-0.0,  A28-0.0; 

float  Bl-0.0,    B2-0.0,   B3-0.0,  B4-0.0,   B5-0.0,   B6-0.0,   D7-0.0,   POO.d. 

float  B9-0.0,   BlO-0.0,  Bll-0.0,  B12-0.0,  B13-0.0,  B14-0.0,  D15-0.0,  niO-0.0; 

float  B17-0.0,  Bie-0.0,  B19-0.0,  B20-0.0,  B21-0.0,  B22-0.0,  P23-n.O,  \\:\    Vi.^.-. 
float  B25-0.0,  B26-0.0; 

float  Cl-0.0,    C2-C.0,   C3-0.0,  C4-0.0,   C5-0.0,   C6-0.0,   C-O.n,   C"  o.o. 
float  C9-0.0,   ClO-0.0,  Cll-0.0; 

float  Dl-0.0,    D2-0.0,   D3-0.0,  D4-0.0,   D5-0.0,   D6-0.0,   D7-0.(i,   n'^-'^.i-, 

float  El-0.0,    E2-C.0,   E3-0.0,  E4-0.0,   E5-0.0,   E6-0.0,   £7-0.0,   f:'^.-0.O; 
float  E9-0.0,   ElO-O.O; 

float  Zl-0.0,    Z2-0.0,   23-0.0,  Z4-0.0,   Fl-0.0,   F2=0.0; 

dydx(4)  -  y[5); 

dydx(6]  -  y(7); 

zeta  -  rho*y  [4]+rhop*y  [6] ;  dydxzeta  -  rho*y  (  5] -•  i  hop' y  ;  7  )  ,- 

Al   -  -(A-C)*y(2]*y|3);  A2   -  J3p*sigma • y ( 2 ] ; 

A3   -  2*MT*zeta*dydx2eta*y [1] ;  A4   -  MT*SQ(zeta); 

A5   -  -MT*SQ(zeta) 'y [2] 'y 13] ;  A6   -  -2*m* (zeta+L2 ) ♦ y [ 4 ) ; 

A7   -  2*m*(2eta  +  L2) 'y (4)*y[2]*y(3);  A8   -  -2*m* (zeta  +  L2)  *y  [  5)  •  >  U ] ; 

A9   -  -2*m*y  (4)  •dyd.x:eta*y  (1)  ;  AlO  -  2  *m*y  I  4  ]  *y  (  5  )  *y  i  1)  ; 

All  -  m»SQ(y(4]);  A12  -  -m*SQ (y ( 4 ] ) 'y | 2 j ' y ! 3 ) ; 

A13  -  -m*yl4)*a;  A14  -  -m*y [ 4) *a»y 1 1 ] *y I  2  J ; 

A15  -  -2»mp* (zeta-Ll) *y (6) ;  A16  -  2*mp* (zeta-Ll ) *y ( 0) 'y 1 2 ) ' y | 3)  ; 

A17  -  -2»mp*  (zeta-Ll)  •y|71*y(l)  ;  A18  -  -2*nip*y  |  6] 'dydxzotcTy  |  1  )  ; 

A19  -  2«mp*y(6] *y[7]*y(l] ;  A20  -  mp*SQ(y|6]); 

A21  -  -np*SQ(y[6))»y(2]«y[3);  A22  -  -mp*y  (  6  )  *ap*cos  (.s  iqT.T  .v)  ; 

A23  -  -rr.p'y  |6]  *ap*cos  (sigma*x)  *y  (1)  *y  12)  ;  A24  -  mp'ap*.';  i  :i  (s  i  rj::!  *  :•: )  ; 

A25  -  mp"'ap*sin{sign'a*x)  *y(6)  *  (SQ(y  (3)+sigina) -SQ(y|2)  )  )  ; 

A26  -  A+A4+A6+A11+A15+A20;  A27  -  A13+A22; 

A28  -  Al-t•A2+A3  +  A5+A7+A8+A9+A10  +  A12+A14+A16+A17+A18  +  A19^A21^A23^/^?^..• 

Bl   -  -{C-A)*y[i;*y(3];  B2   -  -  J3p*  sigm.n  •  y  |  1 )  ; 

B3   -  2*KT*zeta*dydxzet:a*y(2) ;  B4   -  MT*SQ  (zoLa)  ; 
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B5 

B7 

B9 

Bll 

B13 

B15 

B17 

B19 

B21 

B22 

B24 

B26 

CI 

C3 

C5- 

C7 

C9 

CIO 

Dl 
D3 
D5 
D7 

El 
E3 
E5 
E7 
E9 

Zl 
Z2 
Fl 
Z3 

Z4 

F2 


MT*SQ{ 

-2*m*  ( 

-2*m*y 

m*SQ(y 

-m*a; 

-2*mp* 

-2*mp* 

mp*SQ( 

-mp*ap 

-mp*ap 

mp*ap* 

B1+B2+ 


2eta)*y[l)*y(3]; 
zeta+L2) *y|4) •ylll*y|31 
(4) •dydxzeta'y  [2] ; 
I  4 ) )  ; 

(zeta-Ll) *y(6); 

(zeta-Ll)*yI7)*y[21; 

y(6)); 

*cos  (sign)a*x)  ; 


B6  -  -2*m* (zeta»L2 
B8  -  ~2*in*  (zeta*L2 
BIO  -  2*m*yl41*yl51 
B12  -  m*SQ(y|4) )*yl 

B14  m*a*y[4) * (SO 

B16  -  -2*mp* (zeta-L 
B18  -  2*mp*y[6J*y|2 
B20  -  mp*SQ(y|6)) *y 
B23  -  -mp*ap*sin  (si 


*cos(sigina*x)*(SQ(y(3)+sigma)-SQ(ylll  )  )  *y  IGl.- 
slh  (sigrna^x)  *y  (  6)  •y  [1]  *y  [2) ;  B25  -  A« 
B3+B5+B7+B8+B9+B10+B12+B14+B16+B17+B18+B20+B22+ 


)*V 

•yl 
U* 
(yl 
1)  * 
1*  ( 
in 

cirm 

n2A 


ii; 

?); 
yl^ 

3)) 

yp 
•yl 

•y ) 


•  y  I  r  1  ; 

-'■OCyl 
•  y  I  1 1 

\   -HyH:-. 

•  y  I  f; !  ; 
■  1  ]  nn 


1 ) ) )  ; 
*yri 

;'M  .1) 


HI 


-  m*a 


2*m*a*yI51*y(ll; 
yI4]*y(2)»y[3] 


C2 
C4 


-  -m*a*y [4] ; 

-  -2*mp*ap*sin (sigma 


X )  •  V I  <:  ] 


m"a"y  1 1 J 'y  (i:  J 'y  I J  J ;  C4  -  -2*mp*ap*sin  (sigma 

-mp*ap*sin  (sigina*x)  *y  (6) ;  C6  -  -mp*ap*sin  (sigma*x)  *  y 
-2*mp*ap*cos  (sign>a*x)  *y  (7)  *y  [11 ;  C8  -  -mp*ap* cos  (sigma *:< 
mp*ap*cos ( Sigma *x) *yl6]*y(21*y(3); 

C2+CB; 


)  '  y  I  7  )  •  V 


C8  -  -mp*ap* 

1  -  C1+C3+C4+C6+C7+C9; 


yli 
I  f.  1  ; 


M 


CI 


m*(l-rho);  D2   -  -mp*rho; 

-m*a;  D4   -  m*a*y ( 1 ) *y I  3] 

-mMSQ  (y  1 1]  )  +SQ  (y  [2] )  )  •  (y  M 1  *  (1-rho)  -L2-rhop«y  (  6]  )  ; 
k*yl4);  D8   -  D4+D5+D6+D7 ; 


P"" 


'yIM 


-in*rhop;  E2   -  mp*  (1-rhop)  ; 

mp*ap*sin  (sigwa*x)  ;  E4  -  mp*ap*sin  (sigma*x)  *y  |2  ]  *  (y  I  3  )  '  ?*  r-.  i-im  i) 

-mp*ap*cos  (sigina*x)  ;  E6  -  mp*ap*cos  (sigma*x)  *y  [  1  ]  *  (y  !  31  •  2*  r- irjir.i) 
-mp*(SQ(y(l])-fSQ(yI2]))*(y(6)*(l-rhop)+Ll-rho*y(4))  ;  F''  -  rp'yl/l 
kp*y[6];  ElO  -  E4+E6+E7+E0+E9; 

(-(A2  7*B13) /(A2  6*B23)-E1/E3)/((A27*B25)/(A2  6*B2  3)+E5/E3) ; 
(-(A27*B21)/(A26*B23)+A24/A26-E2/E3)/{(A27*B25)  /(A26*B23)  ■•F.5/F.3)  ; 
(-(A27*B2  6)/{A26*B23)+A28/A26-E10/E3)/<  (A27*B25)/  (A26  ^02  3  )  -•  E5/F.3)  ; 
( (C*B13)/(C10*B23)-{A27*B13)/(A26*B23))/ 

( {A27*B25)/ (A26*B23)-(C*B25)/{C10'n2  5)+CS/C10) ; 
( (C*B21)/(C10«B23)-(A27*B21)/(A26*B23)+A24/A2  6)  / 

(  (A27*B25)/(A26*B23)-(C*B25)  /  (CIO '023)  -tCS/CJC)  ; 
( (C*B26)/(C10*323)-C11/C10-(A27*B26)/(A26*B23)+A28/A2  6)/ 

(  (A27*B25)/(A26*B23)-(C*B25)/  (C10*L^7  3)  •C"^/C10)  ; 


dydx(5]  -  ((-r2-D8/D3)/(Z4+D2/D3)+(Fl+D8/D3)/(Z2+D2/D3))/ 

((-Z1-D1/D3)/  (Z2+D2/D3)  +  (Z3-«Dl/ni)/  ( r  •,  i|)?/l>^)  )  ; 

dydx[7)  -  ( (-F2-D8/D3)/(Z3+D1/D3)  +  (F1+D8/D3)/(Z1+D1/D3)  )/ 

((-Z2-D2/D3)/(Zl+Dl/D3)  +  ( 7.1^02/03)  /  (  ".  3n-)l  •  H--.)  )  ; 

dydxI3]  -  (dydx(5]* (C5*B13)/(C*B25)+ 

dydx(7]  •  (  (C10*A24)/  (C*A26)  ■»-{C5*B21)  /  (C*B25)  )  + 

(C10«A28)/(C*A26)+(C5*B26)/(C*B25)-C11/C)/ 

(1-(C10*A27)/(C*A26)-(C5*B23)/(C*B25)); 

dydx[l)  -  dydx(3]*(-A27/A26)+dydx(7]*(-A24/A26)+(-A28/A26) ; 

dydx[2)  -  dydx|3)  *  (-B23/B25) +dydx  (5)  *  (-B13/B25) +dydx  (7]  *  (-B21/B2:<) -P2r,/n;  r.: 
} 


/*•• ***************    MAIN  PROGRAM  *•****•**•******•*«•••• 

main () 

i 

int    i,  j,  nbad,  nok; 

int    11-1000,  12-20C0,  13-3000,  14-4000,  15-5000,  16-6000,  17-700.-; 

int    18-8000,  19-9000,  110-10000; 

float  eps,  hi,  hmin,  dxsavd,  xl,  x2,  xscale,  •ystart; 

float  Hl-0.0,  H2-0.0,  H3-0.0,  Hlp-0.0,  H2p-0.0,  H3p"=0.0,  zcm-0.0,  T:qmn 


/ 
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float  Tl-0.0,  T2-0.0,  T3-0.0,  T4-0.0,  T5-0.0,  T6-0.0,  T?-0.0, 

float  T9-0.0,  TlO-0.0,  Tll-0.0,  T12-0.0,  T13'0.0; 

float  zcmrigid,  It,  Itp,  Ittotal,  Is,  Isp,  Istotal,  wn,  l.nmbdT 

float  •k.e,  *kedel,  *etotal,  •etotaldel,  *h,  *hdel,  "tliota; 

float  "Ec,  *Ecdel,  *Ecp,  *Ecpdel,  *Echype,  'EchypGdol,  •Ecpliyp'-,  * 

float  Q,  Qp,  *eta,  *etap,  Ecmin,  ecfactcr,  ecpfactor,  hsq,  Ecfim! 

float  signec,  signecp,  signq,  signqp,  arg,  argp,  Ecf,  Ecfp; 

float  •Ecexp,  *Ecpexp,  *Ecexpdel,  •Ecpexpdel,  *etali,  *otahp; 


.t:'|.  111-.; 

•  Fcphypfl"  \ 


ystart  -  vectorr (1,N) ; 

ke  -  vectorr (1,MAXARRAY) 

kedel  -  vectorr (1,MAXARRAY) 

etotal  -  vectorr (1,MAXARRAY) 
etotaldel  -  vectorr (l.MAXARRAY) 

h  -  vectorr (l,MAXARRAy) 

theta  -  vectorr (l.MAXARRAY) 

hdel  -  vectorr (l.MAXARRAY) 
Ecexp       -  vectorr (1,MAXARRAY) ; 
Ecpexp      -  vectorr (1,MAXARRAY) ; 


Ec 

Ecpdel 

Echype 

Echypedel 

Ecphype 

Ecphypedel 

eta 

etap 

Ecexpdel 

Ecpexpdel 

etah 

etahp 


vectorr ( 1 ,  "A 
vector  r ( 1 , MA 
vector r (3 ,  MA 
vectorr  ( 3 ,  ''.''■ 
vectorr ( 1 , M" 
vectorr  ( 1 ,  M,". 
vectorr ( i , Mm 
vectorr  ( 1 ,  :','■. 
vectorr ( 
vectorr  ( 1  , .' 
vectorr ( 1 , : 
vectorr ( 1 , ; 


r: 


;•:A^I^■^v) 

XAl'PAV) 
XARl'AV) 
XAPI-AV) 
XAIJPAV) 

xa;u.-ay) 

XAPRAY) 
XAP.PAV) 

■/.^nn^■^i) 

XAI-.PAV) 
XAIU-.AV) 
XARKAV) 


scanf  ("%f%f%f%f%f%f*f%f%f%f%f%f%f%f%f%f%f%C%f%f%r%flf%f'?.f°6f?.f!f -.l", 

4x2,4Il,4l3,tM,Hn,4a,tk,6c,SL,4Ilp,Sl3p,SMp,Smp,Sap,&>-.p,'.  rp, 
4w3p,  Systart(l],tystart(2],tystart(3].iystart('1j,&y2tarv  [5], 
Systartiej.tystartlT],  teps,  tecfactor,  iecpfactor,  4Ecf ,  f.E-rp)  : 

xl-0.0001;       hl-0. 00001;      hmin-1 . Oe-11;    kmax-100000 ; 
mb  -  m;  mbp  -  mp;        12  -  II;         I2p  =  Up; 

Sigma  -  w3p-ystart (3) ; 


/••**  critical  damping  -  input  damping  -  1000.0,  makes  it  crit. 
if  (c  —  1000.0)  c  -  1.5*sqrt  (4*n*k) ; 
if  (cp  ~  1000.0)  cp  -  1.5«sqrt (4*mp*kp) ; 

/•*•  establish  time  Interval  for  saving  data  ***/ 
dxmin    -  l.Oe-4; 
dxlimit  -  (x2-xl) /1750.0; 
dxsav  -  (x2<20.0)  ?  dxmin  :  dxlimit; 


MT  -  M+Mp+4 .0*mb+4 . 0*mbp; 

V  -  (Mp+4 .0*mbp) /MT; 

A  -  Il  +  Ilp  +  2.0*ri)*SQ(a)+2.0*mbp*SQ(ap)  *  (Mp+I  .0*mbp)  *  (1  •  0-v  )  ■ ':  ( ;,)  ; 

C  -  l3  +  I3p+4  .0*rTi)*SQ(a)+4.0*mbp*SQ(ap)  ; 

J3p  -  I3p+4 . 0*mbp*SQ{ap) ; 

rlno  -  m/MT; 

rhop  -  mp/MT; 

LI  -  L* (M+4 .0*mb) /MT; 

L2  -  L*  (Mp+4.0*:T^p) /MT; 

zcmrigid  -  ( (Mp+4 . O'mbp) "L) /MT; 

It  -  Il+M*SQ(zcmrigid)+mb*  (4  .0*SQ(2cmrigid) -'2  .0*SO(a)  )  ; 

Itp  -  Ilp+Mp*SQ(L-zcmrigid) +mbp* (4 . O'SQ (L-zcmr igid)  ^2 . 0* SC  ■-:-)) ; 

Ittotal  -  It+Itp; 

Is  -  13+4 .0»mi*SQ(a) ; 

Isp  -  I3p+4.0*mbp»SQ(ap) ; 

Istotal  -  Is+Isp; 

odeint  (ystart,N,xl,  x2,eps,hl,  hmin,  4nok,  S.-ibad,  derivs,  rkqc)  ; 

for  ( j-1;  j<-kount;  j-t-+) 
1 
zcm  -  {  (Mp+4.0*mbp)  ♦L+m*yp|4]  [  j]-»T,p*yp(6]  (  jl) /MT; 

HI   -  (Il+M*SQ(zcm)+m*  (2.0*SQ(a)+4.0*SQ(zcm)  tSQ(ypM)  I  j;  •■ 
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-2.0*2cn\*yp[4]  I  j]))*yplll  1  j] -m*a*yp(  4  ]  I  j]  *yp  I  3)  [  j )  ; 

H2   -  (I2+M*SQ(2cm)+m* (2.0*SQ(a)+4.0*SQ(zcm)+SQ(yp(4) I j)) 

-2.0*zcm*ypl4] ( j) ) ) *yp|2) | j) -m*a*ypl 6] I j] ; 

H3  -  -m*a*yp[4)  [  j]*yp(l)  (j)  +  {I3+4.0*tn*SQ(a))*yp(3]  I  j); 

Hip  -  (Ilp+Mp*SQ(L-zcm)+inp*(2.0*SQ(ap)+4.0*SQ(L-zcm)+SQ(ypl6]  [  j])^ 
2.0*(L-zcm)*yp(6] ( j)) )*yp(l) ( j) -mp*ap*cos (sigma*xp| j] ) • 
yp(6) I j)*yp[3) [ j]-mp*ap*sigma*cos (sigma*xp| j) ) * 
•  yp(6] [j)+n?>*ap*sin(sigma*xplj])*yp[7) (j); 

H2p  -  (I2p+Mp*SQ(L-2cm)+mp* (2 .0*SQ(ap) +4 . 0*SQ{L-zcm) +50 (yp( 6] I j) ) i 
2.0*(L-zcm)*ypI6] I j]) )*yp(2] ( j] -mp*ap*sin (sigma*xpl j] ) * 
ypl6) ( j]*yp(3] (j]-mp*ap* Sigma* sin (sigma*xpl j] ) * 
yp(6]  ( j]-irp*ap*cos(sigma*xp(  j))  •yp(7]  [j]; 

H3p  -  -mp*ap*cos(sigTna*xp(j])*yp[6)  ( j)*yp(l]  [  j]-mp*ap* 

8in(3igma*xpl  j])«ypI6)  (  j)*yp[2]  (  j)  + (l3p-t-4  .0*mp*SO(ap)  )* 
(ypt3] (jj+sigma); 

h[j]  -  sqrt(SQ(Hl+Hlp)+SQ(H2+H2p)+SQ{H3+H3p)) ; 

argument  -  (H3+H3p) /h ( j) ; 

if (argument>0. 99999999)    argument   -   1.0; 

theta(j]    -   57,2957795131*acos (argument) ; 

hdellj]    -  h[j]-h[l); 

Tl      -   0.5MIl*SQ(yp[l]Ij])+l2*SQ(ypl2)  [j))+I3*SQ(yp[3]  (  j))); 

T2      -   0.5*(Ilp*SQ(yp(l) [jl)+l2p*SQ(ypl2) [J))+I3p* 

SQ(yp(3]  1  j)^r.i'jm->)  ), 

T3      -   0.5*mMSQ(ypl5)  lj))  +  (SQ{yp[l]  lj))+S0(yp|2)  I  j)))*SO(ypMl  I  iD 

-2.0*a*yp(l]  lj)«yp(3]  (jJ'ypMl  I  i)), 

T4  -  0.5*mpMSQ(yp(7]  (  j) ) +SQ  (yp[  6)  ( j) )  MSQ  (ypU)  (  j] ) +SO(yp  12)  I  jl  ) ) 
t-2.0«ap*sin(sigma*xp[  j))*yp[2)  (  j) 'ypl  3J  I  j ) 'yp  I  6  )  U  ] 
-2.0*ap*cos(sigma*xpl jj)*yp(l] [ j) 'ypIS) I j) 'ypl G) ( jl ) : 

T5      -  -m*a*yp[2](j]*yp(5) [j]; 

T6   -  -mp* (-ap*yp[7) ( j] *yp[2] ( j]*cos(sigma*xpl j) ) 
*ap*yp[7) [ j)*ypll) [ j] *sin (sigma*xp 1 j) ) 
-ap*sigma*yp[l] [ j) *yp[6) ( j] *cos (sigma*xp| j] ) 
-ap*sigma*yp(2] [ j) *yp[6] | j) 'sin (sigma'xp ( j] ) ) ; 

T7   -  (0.5*SQ(m*yp[5) ( j)+mp*ypI7] (j)))/MT; 

T8      -    (0.5*{SQ(yp[l] [j))+SQ(yp[2] Ij))) 

*SQ(m*yp(4)  [  j)  ■•mp'yp  (  6]  |  j  )  )  ) /til 

T9      -    0.5*{SQ(ypll] [j])+SQ(ypl2) [ j)))*(M+4.0*m)*SQ{L)* 

SQ( (Mp+1 .0*mp)/nT) , 

TIO   -   0.5*(SQ(yp[l] (j])+SQ(yp[2) I j]))»(Mp+4.0*mp)*SQ(L)* 

SQ( (M+l .0*m) /MT)  , 

Til   -    (m+mp)*(-yp(5] lj])*(m*yp(5] [j]+mp*yp[7] (j])/MT; 

T12   -   -m*yp(4]  [  j)MSQ(yp(l]  Ij))*SQ(yp(2)  (j)))* 

(m*yp[4]  I  j)+mp*ypl6)  [  j)+  (MpM  .0*mp)  •],)  /tlT, 
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T13    -    -mp*yr[6]  ( j]* (SQ(yp(ll  ( j))-SQ(ypl2)  I  j)))* 

(m'ypll)  l3l-tmp*ypl61  I  3l-(f.'1  .0«m)  *  I.)  /Ill; 

ke[j)  -    Tl+T2  +  T3+T4+T5+T6-tT7tT8«T9+T10  +  Tll-tTl?  <ri  i; 

kedellj)  -   ke(j]-ke[l); 

etotal[j]  -   kelj)+0.5*k*SQ(yp[4)  1  j])+0.5*kp*SQ(ypK,l  I  ,1)  ; 

etotaldeKj]  -  etotal  [  j] -etotal  (1 )  ; 

Ec[j}  -   0.5* (Ittotal* (SQ(ypll] ( jl)+S0(ypl2] I jl) ) 

+  lstotal*.<;0(vpni  I  il)  I 

Ecdellj)    -   EC  [j] -EC  ID; 

Ecplj]        -   0.5MIttotalMSQ(ypll]  (  jl)+S0{yp(2]  I  j))) 

+  lstotal*.SO(ypl  ■'1111  •^i<T"i)  ) 

EcpdeKjl    -  Ecp(j]-Ecp|l]; 
) 

wn  -    (Is*ypl3]  (l)+IspMypl3]  ID+sigma)  )/Ittotal; 

lanvbda   -   wn-yp(3](l]; 

lambdap   -   wn- (ypl3)  ( 1  ] +sigTT\a)  ; 

hsq   -    SQ(Ittotal)  MSQ(yp[l]  Il])+SQ(yp(2]  (D)  ) 

+SQ(ls*yp|3)  (U+Isp*  (ypPl  I  1  )  ".-.innn)  ) 

if  (theta [kount] -thetall]  <  0.0)         /....  stable  condition  ••••/ 

{ 

if  (lambda  <  0.0)        (  signec  -  1.0;  signq  "  1.0;  I 
else  {  signec  -  -1.0;  signq  -  -1.0;  ) 

if  (lambdap  <  0.0)       {  signecp  -  1.0;  signqp  -  ).0;  ) 
else  {  signecp  -  -1.0;  signqp  "  -1.0;  I 

) 
else  /.»♦♦  unstable  condition  ♦••'/ 

( 

if  (lambda  <  0.0)        (  signec  -  -1.0;  signq  =  1.0;  ) 

else  {  signec  -  1.0;  signq  -  -1.0;  1 

if  (lambdap  <  0.0)       {  signecp  -  -1.0;  signqp  =  l.O;  ) 
else  (  signecp  -  1.0;  signqp  =  -1.0;  ) 

) 

Ecfinal  -  (Ecf  ~  0.0)  ?  Ec(kount]  :  Ecf; 
Ecpfinal  -  (Ecfp  ~  0.0)  ?  Ecplkount)  :  Ecfp; 

for ( j-1; j<-kount; j++) 
( 
Echypelj]     -  Ec( 1] *Ecf inal/ (Ec (1) + (Ecf inal-Ec ( 1 ) ) * 

exp(ecfactor*Ec£iml*:<p|  j!)  ) 

Echypedel(j]   -  Echype [ j] -Echype ( 1] ; 

Ecphypelj]     -  Ecp [ 1 ] *Ecpf inal/ (Ecp( 1 )  +  (Ecpf inal-Ecp 1 J ] ) * 

exp(ecpf actor  *Ecf  f  im  1  *  :-.p|  j  |  )  ) 

Ecphypedel t j]  -  Ecphype ( j) -Ecphype [ 1]  ; 

Q  -  2.0*Echype( j] Mstotal* (1 .0-(lstotal/lttotal) ) 

-hsq* (Istotal/Ittotal) * (1 .0- (Istotal/Itto' nl) ) 
+ (Istotal/Ittotal) *SQ(Isp) *SQ(sigma) ; 

Qp  -  2.0*Ecphype(  j)*Istotal*  (l.O-dstotal/ILtot.Tl)  ) 

-hsq*  (Istotal/Ittotal)*  (1 .0- (Istotal/ It  Lo'^n  1 )  ) 
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♦  (Istotal/Ittotal)  *SO(Is)  •SO<''iy'"T)  : 

arg  -    ( (Ittotal/ (Ittotal -IsLota) ) ) 

•  (  (Isp*sigma)  ♦  {signq*3'iil  (0)  )  )  )  /  ('  H'  '  (!'•  ■!>  > 
if    (arg    >    0.99999999)    arg-1.0; 

etahlj)  -    57.29S7795131*acos(aig) ; 

argp  -    ( (Ittotal/ (Ittotal-Istotal)  ) 

•  (  (-Is*sigma)  •  (signqp'sqrt  (Op)  )  )  )  /  <•■'{■  '  (•'■    i>  > 
if    (argp   >   0.99999999)    argp-1.0; 

etahplj)  -    57.2957795131*acos(argp) ; 


Ecexp[j]  -       (Ec{l)-Ecfinal)  •exp(ecfactor*xp(  j) )  ^Fcf  iml  . 

EcexpdeKj]        -   Ecexp(  j] -Ecexp(  1)  ; 

Ecpexpl  j]  -    (Ecpll)-Ecpf  inal)  *exp(ecpfactor  *xp|  j  J )  iF.'-pf  i  n-i  1  ; 

Ecpexpdel(j;    -  Ecpexp( j) -Ecpexpl 1] ; 

Q  -   2.0*Ecexp(  j)  •Istotal*  (l.O-dstotal/Ittoi  Tl )  ) 

-hsq*  (I stotal/ Ittotal)  •  (1 .0- (Istotal/ Itt.ot.n  1  )  ) 
♦  (Istotal/Ittotal)  *SQ(Isp)  ♦SQ(sigm;:»)  ; 

Qp  -   2.0*Ecpexp(  j]*Istotal*  (l.O-dstotal/UtoL.Tl  )  ) 

-hsq*  (Istotal/Ittotal)  *(1.0-(Istotal/Ut<-i.-il)) 
+  (Istotal/Ittotal)  "SOds)  *SQ(sigma)  ; 

arg  -  ( (Ittotal/ (Ittotal-Istotal) ) 

•  (  (Isp*sigma)  +  (signq*sqrt  (0)  )  )  )  /  ('-T  •  (li'-'i)  ) 
if  (arg  >  0.99999999)  arg-1.0; 

etalj]       -  57.2957795131*acos(arg); 

argp        -  ( (Ittotal/ (Ittotal-Istotal) ) 

*  (  (-Is*sigma)  +  (signqp*sqrt  (Qp)  )  )  )  /  (r.qrL  (li-'i)  ) 
if  (argp  >  0.99999999)  argp-1.0; 

etaplj)      -  57.2957795131*acos (argp) ; 
) 


initialize  graphics  **•*«•*•*••/ 


init  0 ;   color_scale  "cyanblue") ; 

grey_scale("greyscale:") ;   windowO () ; 

bgcol(7);   erase  0;   :olor(0); 

move (220, -110);  print: ("KINGORI  DUAL  SPINNER  MASS  SPRING  SYSTLn"); 

windou(22, -85, 979,872  ;   bgcol(3); 

erase  0;   scale (0, 10*::.  10000, 0) ;    rect (0, 0, 10000, 10000) ; 

sized,  1)  ; 

rrve  (20, 110  +  200)  ;    dateO; 

rrve (-160, 110+70) ;   printf ("0") ; 

rove (10, 110-90) ;     printf ("w3  -  w3  initial";; 

vector(0, IS. 110,19) ;  rrve (-160, 19+70) ;   printf ("0") ; 

-ove (10, 19+200) ;  printf ("w2") ; 

vector (0,13,110.18) ;  r-ve (-160, 18  +  70)  ;   printf ("0") ; 

-ove (10, 18+200) ;  printf ("wl") ; 
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vector  (0,  1"?. 110. 17) 


move  (-160, 17*70) 
movedO,  17  +  200) 


printf ("0") ; 
prinLfC'z  prime  <<o(.") 


vector (0, 16, 110,16) ;  move (-160, 16+70) ;  print f ("O") ; 

move (10, 16+200) ;  printf ("2  prime"); 

vector (0, 15,110,15) ;  move (-160, 15+70) ;  printf ("0") ; 

move (10, 15+200) ;  printf ("z  dot"); 


vector (0, 11,110,14) 


vector (0, 13, 110,13) 


move (-160, 14+70) ;   printf ("0") ; 
move (10, 14  +  200) ;  printf ("z")  ; 

move (-160, 13+70) ;   printf ("0"); 
move  (10, 13  +  200)  ;  printf  (")ce-ke  i  (r), 


o  -  f"  i  f  I  > )  " ) ; 


vector (0, 12,110,12) ;  move (-160, 12+70) ;   printf ("0"); 

move (10, 12+200) ;  printf ("Ec-Ec  i  (r),   Echypc-Echypc  i  (b)"); 

vector  (0, 11, 110.11) ;  move (-160, 11  +  70) ;   printf("0"); 

move  (10, 11  +  200) ;  printf ("Ecp-Ecp  i  (r),  Ecphype-Ecphypc  i  (!)"); 

vector (0, 0, 110, 0) ;    move (-160, 0+70) ;   printf("0"); 

move  (10,0  +  200) ;  printf  ("theta,  eta  (b) ,  etap  (r),  h  -  »i  imii.»l") 


vector (0, 
vector (0, 
vector (0, 
vector (0, 
vector  (0, 
vector  (0, 
vector  (0, 
vector  (0, 


19+-750,11, 
19+-500,11, 
19+250,11, 
19-250,11, 
19-500,11, 
18+-250,11, 
18-250,11, 
18-500,11, 


19+750) 
19+500) 
19+250) 
19-250) 
19-500) 
18+250) 
18-250) 
18-500) 


xscale  -  10000 .0/xpIkount) 


vector (19,19+750, 110, 19+750) 
vector  (19, 19+500, 110, 19+-500) 
vector (19, 19+250, 110, 19+2  50) 
vector(19,19-250, 110, 19-250) 
vector (19, 19-500, 110, 19-500) 
vector(19.18+250, 110, 18+250) 
vector (19, 18-250, 110, 18-250) 
vector (19, 18-500, 110. 18-500) 


color  (1) . 
for  (j-1. 


jOcount;  J+-  +  ) 


blue  -  total  energy,  rotor  hype 


( 

vector ( (int) ( xscale •xptj) ) , (int) (13+(0.01*etotaldel ( j)) ). 

(int)  (xscale'xpl  J+-1)) ,  (int)  (13+(0  .01*etotaldel  (  jU)  )  )  ) 
vector  (  (int) (xscale*xp[ j) ) ,  (int)  (13+ (0. l*etotaldel ( jl ) ) , 

(int)  (xscale*xp[ j  +  1]),  (int)  (13+ (0 .l*etotaldel  |  j  +  l) ) ) ) ; 

vector (  (int)  (xscale 'xpl j] ) ,  (int)  (12+(50 .0*Echypedel | jl) ) , 

(int)  (xscale«xpl  J+-11) ,  (int)  (12+ (50.0*Echypedcl  1  j+l)  )  )  ) 

vector  (  (int)  (xscale*xp( j) ) ,  (int) (0+ (100 .0*etah| j) ) ) . 

(int)  (xscale*xp[ j  +  lj) .  (int)  (0+ (100 . 0»etah ( j  +  l ) ) ) ) ; 
vector ( (int) (xscale 'xp ( j] ) , (int) (11+ (1000 .0*etah [ j) ) ) , 

(int)  (xscale-xpl  j  +  l]  )  ,  (int)  (11+ (1000  .  O'etali  (  j  +  l ) )  )  )  ; 

} 


color (5);        /♦«...  purple  -  platform  exponential*****/ 
for  ( j-1; jOcount; j++) 

{ 

vector (  (int)  (xscale*xp[ j) ) ,  (int) (11+ (50.0*Ecpexpdel ( j] ) )  , 

(int)  (xscale'xpl  j  +  l) ) ,  (int)  (11+ (50.0*Ecpexpd'?l  I  j<  )  j 

vector  (  (int)  ( xscale 'xpf j] ),  (int) (0+ (100 .0*etap| j] ) ) , 

(int)  ( xscale * xp I j  +  l) ) ,  (int)  (0+ (100. 0*etnp| jU 1 ) ) ) ; 

vector (  (int)  (xscale "xpl j] ) .  (int) (11+ (1000 .0 •etap ( j) ) ) , 

(int)  (xscale*xp(  J+-1]),  (int)  (11+ (1000.0*  etapl  j  •  1  ))  ) ) 


)  ) 
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) 

color(6);  /•«..•    yellow    -    rotor    exponential*'***/ 

for    (  j-1;  j<)count;  j*-? 

( 

vector  (  (int)  :xscale*xp(  j) ),  (int)  (12^  (50  .0*Ecexp<1o  1  I  jl  )  )  . 

(int;  (xscale*xp(  j  +  11)  ,  (int)  (12+ (50  .  0 'F.ccxpdcl  1  j'  I  1  i  m 

vector ( (int)  ;xscale*xp( j) ),  (int)  (0+ (100 .0*eta  1  j  ) )  )  , 

(inti  (xscale*xp[  j  +  D) .  (int)  (0+ (100  .  O'et.i  I  jO  1  )  )  )  ; 

vector  ( (int.-  ixscale*xp(  j)  ) ,  (int)  ( 11*  (1000  .0*cta  (  jl )  )  , 

(ir.f  (xscale*xpl  j4l])  ,  (int)  (11-t  (1000.0'ctAl  j  i  1 1  )  )  I  : 

1 

color(O);-       /•♦...  blacJc  -  platform  hyperbolic  *•**•/ 
for  ( j-1;  jOcount;  j*-*; 


( 

vector  ( (int)  (xscale*xp(  j)  ),  (int)  (11+ (50  .0*Ecpliypcd'7l  I  jl )  )  . 
(int:  (xscale*xpl  j  +  11)  ,  (int)  (11+  (50  .  0*Ecpliypedel  |  j«  1 

vector  ( (int) (xscale*xp( j] ),  (int)  (0+(100.0*etahplj))), 

(int;  (X3cale*xpl j+lj)  ,  (int)  (0+ (100 . 0*etahpl jO ) ) ) ) ; 

vector  (  (int)  (xscale*xp(  j) ),  (int)  (11  + (1000  .0*etal>p  |  j) )  )  , 

(inti  (x3cale*xp( j+1)) , (int) (11+ (1000 . O'ctohpl j' 1 ) ) ) 

) 


1)  )  )  ; 


color  (4);        /*«»-.  ^g^j  _  rotor  **•*•/ 
for  ( j-1; jOtount; j+*) 
{ 
vector ((int)  (xscale*xp( j) ) ,  (int)  (110+ (14  32 .3914  8763* (yp(3)  ( jl-yp 
(int)  (xscale*xpli  +  l)),  (int)  (110+ (1432 . 39448783* (yp  |3)  [ j  +  ll-ypi 

vector ( (int) (xscale*xp( j) ). (int) (19+ (1432 .394 48 78  3* yp ( 2 ] I j I 
(int; (xscale*xp| j  +  11)  ,  (int)  (19+ (1432 . 394  4  0783*yp| ? ) 

vector  (  (int) (xscale*xp[ jl ) ,  (int)  (18+ (14  32 .394  48783*yp| 1 1  I j) 
(int)  (xscale*xp( j+11)  ,  (int)  (18+ (1432 . 394  4  87  83*ypn 1 

vector  (  (int)  <>:scale*xp[  jl  ),  (int)  (17+ (1000*yp(71  [j])) , 

(int)  ;xscale*xp(  j  +  11)  ,  (int)  (17+ (lOOO'yp  (7  1  (  j-tl]  )  )  )  ; 

vector  (  (int)  i::scale*xp[  jl ) ,  (int)  (16+ (1000*yp|  61  1  j]  )  ). 

(int)  !xscale*xp[ j  +  11) ,  (int)  (16+ (1000*yp | 6 ]  1  j  +  l] ) ) ) ; 

v-  -z{  (int)  .'xscale*xp{  jl).  (int)  (15+ (1000*yp(5}  I  j]  )  ), 

(int)  ;xscale*xp( j  +  11),  (int)  (15+ (lOOO'yp ( 5] ( j«ll  ) ) ) ; 

vector  (  (int)  .>:scale*xp(  jl  )  ,  (int)  (14+ (1000*yp(  4  1  [  j]  )  )  , 

(int;  :xscale*xp( j  +  11) ,  (int)  (14+ (1000*yp M ]  (  j  "  1  )  )  )  )  ; 

vector  (  (int)  (>:scale«xp(  jl ) ,  (int)  (13+ (0  .  01*)ccdcl  (  jl )  ), 

(int)  .xscale*xp( j  +  11) ,  (int)  (13+ (0 .01*kcdel  I  jO ) ) ) ) ; 

vector  (  (int)  (:;scale*xp(  jl ) ,  (int)  (13+ (0  .  l*)(edel  [  jl  )  )  , 

(int)  ::<scale*xp(  j  +  11),  (int)  (13+ (0  .  l*)cedel  (  j-" )  )  )  )  )  . 

vector  (  (int)  (>:scale*xp(  jl  ),  (int)  (12+ (50  .O'Ecdel  I  j  1 )  ). 

(int)  .xscale*xpl j  +  11) ,  (int)  (12+ (50 . O'Ecdel ( j • 1 ) ) ) ) ; 

vector  (  (int) (x3cale*xp( jl ) ,  (int) (11+ (50 .0*Ecpdcl I jl ) ) , 

(int)  :xscale*xp( j+11), (int)  (11+ (50.0*Ecpdel  I  j'lD ) ) 

vector  (  (int)  (>:scale*xpl  jl )  ,  (int)  (0+ (1000  .  0*hdel  (  jl )  )  , 

(int)    xscale'xpl  j  +  11 )  ,  (int)  (0+ (1000  .  O'lidc]  I  j'  1  1  )  )  )  : 


!  -^  1  I  1  1  )  )  ) 
3)  il)  )  )  )  ) 

)  )  , 

Ij  ■  1 i ) ) ) ; 


)  )  , 


» 1  !  )  )  )  ; 
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vector  (  (int)  (xscale*xp|  j  )  )  ,  (  int )  (OM  100  .  0  •Lliot  n  |  j  )  )  )  , 

(int)  (xscale*xp(  j-"!  ))  ,  (int)  (0<  ( 100  .  0*  t  lioL.i  1  j  •  1  1  )  ))  : 

vector  (  (ir.t)  (xscale*xp[  j]  ),  (int)  (11-'  ( 1000.0 *thctn  |  il  )  )  , 

(int)  (xscale*xp|  j  +  D).  (int)  (IH  (1000 .  OM  hot.il  j  i  1  I  )  )  )  : 

vector  (  (int)  (xscale*xp(  j)  )  ,  (int)  (OMO  .0001*ke  (  j) )  ) . 

(int)  (xscale'xpl  j  +  D)  ,  (int)  (0+ (0  .  0001  *kc  1  j<  1  )  )  )  )  . 

vector ( (int)  (xscale*xpl j) ) ,  (int)  (0+ (0 . 0001 'etotal  I  j) ) ) , 

•   (ir.t)  (xscale*xpl  j  +  lj),  (int)  (0*  (0  .  OOOl'ctotnl  (  j '  I  1  )  )  )  ; 

vector  (  (ir.t)  (xscale*xpl  j] ) ,  (int)  (0+ (0  .  l*h  (  j) ) )  - 

(ir.t)  (xscale*xpl  j  +  1])  ,  (int)  (0*  (0  .  l*h  (  j  +  l )  )  )  )  ; 


windowQ  0 ; 
color  (0); 

/••♦   print  initial,  final  data   ' 
move (350, 675) ;  printf {"         INITIAL  CONDITIONS 


koinl 


1.(1",  I'Miiit  ) 


/• 

move (22, 910) ; 

move (500, 910) ; 


printfC'time  -%3.0f  seconds",  xp  (Jcount  J ) 
printf ("eps  -  %e    \d/\d" , eps, nbad, nok) 


move (22, 94  0) ; 

prlntf("Wl  -»4.2f  W2  -\4.2f  W3  -*4 . 2f ", yp ( 1 ) 1 1 ] , yp I  2 )  | 1) , yp I  1 1  ( 1 1 ) 

move  (500,  94  0)  ;  pri.-.tf  (''w3p-  *9  .  7f        L  -  %3  .  If ",  w3p,  L)  . 

move(22,980);  prir.tf{"Il  -*4.0f  12-%4.0f   13  -14  .  Of  ",  1 1 ,  T?,  I  ■«)  ; 

move(22,1010) ;  printf ("M  -  %4.0f  m  -%7 . 5f ", M, m) ; 

move (22, 1040) ;  printf("a  -%4.1f  k  -%4.1f    c  -%4 . If ", a, k, c) ; 

move(22,  1070)  ;  printf("2  -\A  .2£  2dot-\4  .  2f  ",  yp  (4 )  ( 1 ) ,  ypl  .S)  ( 1 J )  : 

move  (500,  980)  ;  printf  ("Ilp-*4  .  Of  I2p-%'1.0f    I3p-\4  .  Of  ",  Dp,  I  7p,  I  :'.p)  ; 
move (500, 1010) ;  printf ("Mp  -*4.0f  mp-%7 . 5f ",Mp,mp) ; 
move  (500,  1040)  ;  printf  ("ap-\4. If  )cp-%4.1f     cp-t4  .  If  ",  ap,  kp,  cp)  ; 
move (500, 1070);  printfC'zp  -%4.2f  zdotp  -%4 .2f ", yp|6) ( 1 ) . ypl 7) | 1 )) ; 


move (22, 1110)  ; 
move (22, 1135)  ; 
move (22, 1160) ; 
move(250, 1160) , 
move(500, 1160) , 

printf  ("e  celper-*5  .2f  ",  100.  0*(  etotal  Ikount] -etotal  1 1 ) )  /'^loI-.t  J  I  1  ) 
move(730, 1160) ;  printf ("Ecp  delper-%5 . 2f ", 100.0* (Ecp[koiint I -Erp ( 1 ) ) /Frp\ 1 1 ) 


printfC'h  i-*9.4f   ke  i-%7.  If  ",  h  ( 1] ,  kel  1  ]  )  ; 
printf ("h  f-%9.4f   ke  f-%7 . If ", h (kount] , kc ( kount )) ; 
printf ("hdelper-*5.2f", 100.0* (h[kount)-h|l)  ) /hi  1)  )  ; 
printf  ("kedelper-\5.2f",  100.0*  (ke  (kount  1 -ke  1 1  1  )  /(-."i  1  )  ) 


move(500, 1110) ;  printf("Ec  i-%9.4f 
.T.ove  (500,  1135);  printf  ("Ec  f-%9.4f 


theta    i-%7.5f",EcIl]  ,  thet-Tl)  1  )  ; 

theta    f-%7.5f  ",  Ecikount  J ,  tlicL.T  I  ^oniii  )) 


move (22, 900)  ; 

printf ( 

" 

pla 

tf  orm 

move (22, 920)  ; 

printf ( 

"Up 

- 

*6 

2f" 

,  Up) 

move (240, 920) ; 

printf ( 

"11 

- 

%6 

2f" 

,11); 

move (22, 940) ; 

printf ( 

"I3p 

- 

*6 

2f" 

,  I3p) 

move (240, 940)  ; 

printf ( 

"13 

m 

i6 

2f" 

,13); 

move (22, 960)  ; 

printf ( 

"Mp 

- 

*6 

2f" 

,Mp); 

move (240, 960)  ; 

printf ( 

"M 

- 

16 

2f" 

,M)  ; 

move (22, 980); 

printf ( 

"mp 

- 

\€ 

2f" 

,mp)  ; 

move (240, 980)  ; 

printf ( 

"m 

- 

%6 

2f" 

,m); 

.-ove(22,1000)  ; 

printf ( 

"ap 

- 

«6 

2f" 

,ap); 

.-ove(240,  1000)  ; 

printf ( 

"a 

- 

%S 

2f" 

,a)  ; 

move (22, 1020) ; 

printf ( 

•kp 

- 

%6 

2f" 

,  kp)  ; 

move (240, 1020) ; 

printf ( 

"k 

- 

%6 

2f" 

.k); 

rotoi ") ; 
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move (22, 1040) ;   printf 
move(240, 1040) ;  printf 


move 
move 

move 
move 
move 
move 
move 
move 
move 
move 


move 
move 
move 
move 
move 
move 
move 
move 
move 
move 
move 
move 


move 
move 
move 
move 
move 
move 
move 
move 
move 
move 

move 
move 

move 
move 
move 
move 
move 
move 


22, 1060) ; 
240, 1060)  , 

22, 1100) ; 
240, 1100), 
22, 1120) ; 
240, 1120) 
22,1140) ; 
240, 1140), 
22,1160) ; 
240,1160), 


printf 
printf 

printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 


480, 
730, 
480, 
730, 
460, 
730, 
480, 
730, 
480, 
730, 
480, 
730, 
'P 


900) 
900) 
920) 
920) 
940) 
940) 
960) 
960) 
980) 
980) 
1000 
1000 
rintf ("edelp 


printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 


"cp 
"c 


\0.2f*,cp) 
\  6  .  2  f  •• ,  c )  ; 


480,1030) 
730,1030) 
480, 1050) 
730,1050) 
480,1070) 
730,1070) 
480, 1090) 
730,1090) 
480,1110) 
730,1110) 

480,1130) 
730,1130), 

480, 1155), 
620, 1155). 
760.1155). 
480, 1175) i 
620,1175); 
760, 1175), 


printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 
printf 

printf 
printf 

printf 
printf 
printf 
printf 
printf 
printf 


"L  -  \6.'IC",  I,)  ; 

"Is/It  -  %6.3C", Istotal/IttoLnl) 

"zp  -  \6.2f",yp|61  II  J) 

"z  -  %6.2f",ypM)  ID) 

"dip  -  %6.2f",yp|7)  (D) 

"dz  -  %6.2C",ypl5)  ID) 

••w3p  -  %6.2C",ypl3J  ll)-»sigm.n) 

••w3  -  %6.2f",ypl3)  (D) 

"wl  -  %6.2f ",yp|l)  HI) 

••w2  -  %6.2f",yp|2)  11}) 


"\3.0f    seconds",  xptJcount ))  ; 

••hi  -    %9.4f,h[l))  ; 

••%d    /    \d  /    %d^^,  nbad,  nok,  kount)  ; 

••h    f  -    %9.4f  ■•,h|kount)  )  ; 

••eps  -    \7e'^,  eps)  ; 

••hdelp  -    \5.3f,  100.0*  (li(kounLl-hn  l)/i-.j  ' 

••ke    i  -   %9.4f,kell))  ; 

••e    i  -   *9.4f,etotalll))  ; 

•'ke    f  -    \9.4f  ",kelkountJ)  ; 

••e    f  -    *9.4f,etotalIkount  J)  ; 

"kedelp  -    %5  .  3f,  100.0*  (ke  [kount] -ko  | )]) /v 

-   %5.3f",  100.0*  (etotal(kount)-etonf>l  ID)  /- 


^\\]) 


-«1  I  i  1) 


"Ecpi    -  \9.3f",Ecp(ll); 

"EC  i    -  %9.3f".EcIl)); 

"Ecpf    -  %9.3f",Ecp(kount]); 

"Eg  f    -  %9.3f",Ec{kountl) ; 

"lambdap  -  %5 . 3f ", lambdap) ; 

•'lambda  -  %5.  3f  ",  lambda)  ; 

"slgnqp   -  %5. 2f ", signqp) ; 

"signq   -  %5 . 2f ", signq) ; 

"dEcp/lambdap-  *6.3f",  (Ecp  (kount ) -Ecp(  1  ) )  /  l.-«T.b'lTp) 

"dEc/Iambda   -  %6.  3f  ",  (Ec  [kount ) -Ec  1 1 )  ) /I iiri^'i-i)  ; 

"Ecfp    -  %8.6f".Ecfp) ; 
"Ecf     -  l8.6f",Ecf) ; 

"theta  i  -%6  .2£", theta 1 1 ) ) ; 
"etap  i  -\6.2f ", etapll) ) ; 
"eta  i  -%6.2f",etall)) ; 
"theta  f  -%6 .2f", theta Ikount )) ; 
"etap  f  -%6.2f ",etap|kount) ) ; 
"eta  f  -%6.2f  ■•,etaIkount)  )  ; 


/••*   upper  right  hand  corner   ***/ 
move(600,-33) ;  printf ("Ec  initial/final  ~   \g   /   %g", Ec ( 1 ) , Ec [kount ] ) : 
move{600,-13)  ;  printf("Ecp  initial/final  -  %g  /  lig",  Ecpl  1  ] ,  Ecp|kc\jr-.:. ;  ) 
move(600, 27)  ;   printf ("Ec  factor  /  Ecf   -  %g  /  %g", ecfactor, Ecf ! ; 
move{600, 47)  ;   printf ("Ecp  factor  /  Ecfp  -   \q  /   *g", ecpfactor , Ecfp) : 


move(9e0, -33) 
move(980,-13) 
move (980, 7) ; 


printf ("%d", nbad) ; 

printf  (•'%d",nok)  ; 

printf ("%d",  kount) ; 


move (980, 27) ;  printf ("nbad") 

move (980, 47) ;  printf ("nok") ; 

move (980, 67) ;  printf ("kount") ; 

move (980, 97)  ;  printf ("%6 . 4 f",  (100.0'  (ke (kount ) -ko ( i ) ) ) /koM 1 ) 
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move  (980,  IIT);   pr  intf  ("»6.  If  ",  (100.0Mct.ol.il  (kount  ) —tol.i  1  1  H  ))/•"<'''  I  H  ) 
move  (980,  137)  ;  pr  intf  ("iG  .4  f  ",  ( 100  .  0*  (l>  (  konnt.  | -I.  I  1  ))  )  /hH  ))  ; 


move (980, 157) ;  printf 
move (980, 177) ;  printf 
move (980, 197) ;  printf 

move (980,  227) ;  printf 
move  (980,  247)  ;  prir.tf 

move(980,267) ;  printf 
move (980, 287)  ;  printf 

move (980, 317) ;  printf 
move  (980,  337)  ;  prir.tf 

move (980, 357) ;  printf 
move (980, 377) ;  printf 

move (980, 407) ;  printf 
move (980, 427) ;  printf 

move (980. 447)  ;  printf 
move (980, 467)  ;  printf 

move (980, 497)  ;  printf 
move (980, 517) ;  printf 
move (980, 537) ;  printf 

move (980, 557) ;  printf 
move (980. 577)  ;  printf 
move (980, 597) ;  printf 

move (980, 627)  ;  printf 
move (980, 647)  ;  printf 
move (980, 667)  ;  printf 

move (980, 687) ;  printf 
move (980, 707)  ;  printf 
move (980, 727)  ;  printf 

move (980. 757) ;  printf 
move(980,  777) ;  printf 


"kedel") 
"edelp") 
"hdelp") 


•*\g",theta(l))  ; 
"\g", thetalkount ) ) ; 

"th  1") ; 
"th  f") ; 

"*g".etall)) ; 
"%g", eta (kountl ) ; 

"eta  1"); 
"eta  f"); 

"%g",etap(l)) ; 
"%g", etaplkount) ) ; 

"etapl") ; 
"etapf ") ; 

••%g",Ec[kountJ-Ec(lJ)  ; 
"%g", lambda) ; 
"%g", signq) ; 

"Ecdel") ; 
"lambda") ; 
"signq") ; 

"%g",Ecp(kount]-EcpIll) 
"%g", lambdap) ; 
"%g", signqp) ; 

"Ecpdel"); 
"lambdap") ; 
"signqp") ; 

"%g",    Istotal/Ittotal) ; 
"Is/It"); 


free 

_ma  t  r  i  X 

free" 

_vector 

free' 

"vector 

free' 

vector 

free" 

_vector 

free" 

_vector 

free' 

vector 

free" 

vector 

free" 

vector 

free" 

vector 

free" 

"vector 

free" 

vector 

free" 

vector 

1 

( yp, 1 , N. 1 , KAXARRAY ) ; 
(xp, 1,  KAXARRAY) ; 
(ystart, 1,N) ; 
(ke, l,MAXARRAy) ; 
(kedel,!,. MAX ARRAY) ; 
(etOtal,l,MAXARRAY) ; 
(etotald€l,l,MAXARRAY) 
(h,  l.MAXAJ^RAY)  ; 
(theta,  1,.»^AXARRAY)  ; 
(hdel,l,.'<AXARRAY)  ; 
(Ecexp,l..*'JiXARRAY)  ; 
(Ecpexp,  1,.«AXARRAY)  ; 
(etah,  1 ,  H-.XARRAY ) ; 


free_ 

vector 

free 

vector 

free' 

vector 

free' 

vector 

free' 

vector 

free^ 

vector 

f  ree_ 

vector 

free 

vector 

free' 

vector 

free' 

vector 

free" 

vector 

free 

vector 

free" 

"vector 

(Ec,  1,  riA.XAPPA-.')  ; 
(Ecdcl,  l.KA.YAi-r-Y)  ; 
(Ecp.  l.HAXAP.PAV)  ; 
(Ecpdel ,  1 ,  KAy.A.f  P.AY  )  . 
(Echype,  1,  t•1AXA^^^Ay  )  ; 

(Echypcdcl ,  1 ,  :iA:-:ArrAY )  . 

(Ecphypc,  1,MAXAP.I'AV)  ; 

(Ecphypcd«l,  1 ,  ;:A:<Ai?i'.Ay ) 

(eta, 1,MAXAPRAV) . 
(etap,  l,MAXAIU'.Ar)  ; 
(Ecexpdel  ,  1,  MAXAPI'AV)  ; 
(Ecpexpdol,  i.  tlAXAKi'.AY)  . 
(Ctalip,  l,MA/.Al'.i-AY)  ; 
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